Solving the heat equation

heat.c

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double init_fluid(double x) {
  return 10.0 * sin(2.0 * M_PI * x / 30.0) + 8.0 * cos(x / 2.0) +
         2.0 * sin(x / 2.0);
}
int main() {
  double u[500], v[500], t = 0.0, k = 50.0, h = 12, delta_t = 0.001;
  for (int i = 0; i < 500; i++) {
    u[i] = init_fluid(((double)i * 30.0 / 500.0));
    v[i] = 0.0;
  }

  for (int frame = 0; frame < 200; frame++) {
    FILE *fp = fopen("pl.dat", "w");
    assert(fp != NULL);
    for (int i = 0; i < 500; i++)
      fprintf(fp, "%f %f\n", (double)i * 30.0 / 500.0, u[i]);
    fclose(fp);
    char str[300] = {0};
    sprintf(str,
            "gnuplot -e \"set term png; plot [0:29] [-15:15] 'pl.dat' with "
            "lines\" > frame-%03d.png",
            frame);
    puts(str);
    system(str);

    for (int ctr = 0; ctr < (int)(0.04 / delta_t); ctr++) {
      for (int i = 0; i < 500; i++)
        v[i] = k *
               (u[i + 1 > 499 ? 499 : i + 1] + u[i - 1 < 0 ? 0 : i - 1] -
                4.0 * u[i]) /
               pow(h, 2.0);                                                                   for (int i = 0; i < 500; i++)
        u[i] += delta_t * v[i];
    }
  }

  return 0;
}