// Solve 2nd order differential equation using Euler method
// C
// compile as ::   gcc -o DE2nd.exe DE2nd.c -lm 
// execute as ::   ./DE2nd.exe
// output files :: output.analytic     t  x  
//                 output.numerical    t  x  v
#include <stdio.h>
#include <math.h>
//
#define X0 3.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 dxdt){  
    double dx2dt2 ; 
    dx2dt2 =  -2.0 * dxdt - 10.0 * x;
    dx2dt2 =  -4.0 * x;
    return dx2dt2;
}
// ----------------------------------------------------------
// analytic solution 
double sol(double t){  
    double solution ;
    solution = X0 * exp(-t) * cos(3.0 * t);
    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,xnew,vnew;
  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
      xnew = x    + dt * dxdt ;
      vnew = dxdt + dt * dvdt(x, dxdt);
      x    = xnew;
      dxdt = vnew;
    // next t
      t += dt;
  } // end of t-loop
  // close files 
  fclose(fp1);
  fclose(fp2);
}