#include <stdio.h>

int main() {
  double t, x, y, w2, A;
  double dx, dy, dt;
  int i, n;
  FILE *dat, *gp;
  
  A = 1;
  w2 = 1;
  dt = 0.5;
  n = 40;
  dat = fopen("HarmonicOscillator.dat", "w");

  // initial conditions
  x = A;
  y = 0;
  t = 0;
  for (i=0; i< n; i++) {
    dx = y * dt;
    x += dx;
    dy = -w2*x*dt;
    y += dy;
    t += dt;
    fprintf(dat, "%f %f %f\n", t, x, y);
  }
  fclose(dat);
  gp = popen("gnuplot", "w");
  fprintf(gp, "plot 'HarmonicOscillator.dat' u 2:3 w l\n");
  fflush(gp);
  getchar();
  fclose(gp);
  return(0);
}