// Solve 2nd order differential equation using Euler method
// C
// compile as ::   gcc -lm -o DE2nd.exe DE2nd.c
// execute as ::   ./DE2nd.exe
// output files :: output.analytic     t  x  
//                 output.numerical    t  x  v
#include <stdio.h>
#include <math.h>
//
#define X0 1.0       /* Initial Value  x(0) */
#define V0 0.0       /* Initial Value  v(0) */
#define T0 0.0       /* starting time t0 */
#define TMAX 10.0    /* ending time tmax */
#define OUTPUTSTEP 10 /* output data every xx steps */
// ----------------------------------------------------------
// equation of motion
//     dvdt = right-hand side of the 2nd order DE
double dvdt(double x){  
    double f =  -4.0 * x;
    return f;
}
// ----------------------------------------------------------
// analytic solution 
double sol(double t){  
    double solution = X0 * cos(2.0 * t);
    return solution;
}
// ----------------------------------------------------------
int main(void)
{
  char filename1[] = "output.numerical";
  char filename2[] = "output.analytic";
  FILE *fp1, *fp2;
  double dt=0.01; // delta t
  double t,x,dxdt;
  int icount=0;
  // open files
  fp1 = fopen(filename1, "w");
  fp2 = fopen(filename2, "w");
  // initial set up
  t = T0;
  x = X0;
  dxdt = V0;
    printf("dt= %8.4f \n",dt);
    printf("     t       numerical     analytic     diff \n");
  // t-loop
  while(t < TMAX){
    // check accuracy and output
      if(icount % OUTPUTSTEP == 0){
          // output
          printf("%10.3f %11.5f %11.5f %12.8f \n", t,x,sol(t),x-sol(t));
          fprintf(fp1,"%12.5f %12.5f %12.5f\n", t,x,dxdt);
          fprintf(fp2,"%12.5f %12.5f\n", t,sol(t));
      }
      icount += 1;
    // Forward Difference
    dxdt += dt * dvdt(x);
    x    += dt * dxdt   ;
    // next t
    t += dt;
  } // end of t-loop
  // close files 
  fclose(fp1);
  fclose(fp2);
}