// Solve 1st order differential equation using Euler method
// C
// compile as ::   gcc -lm -o DE1st.exe DE1st.c
// execute as ::   ./DE1st.exe
// output files :: output.analytic     t  x  
//                 output.numerical    t  x
#include <stdio.h>
#include <math.h>
//
#define X0 2.0       /* Initial Value  x(0) */
#define T0 0.0       /* starting time t0 */
#define TMAX 10.0    /* ending time tmax */
#define OUTPUTSTEP 10 /* output data every xx steps */
// ----------------------------------------------------------
// differential equation
//     dxdt = right-hand side of the 1st order DE
double rhs(double x, double t){
    double dxdt ;
    dxdt =  -0.5 * x + exp(-0.2 * t) ;
    dxdt =  -0.5 * x + sin(t) ;
    dxdt =  -0.2 * x * t;
    dxdt =  -0.5 * x ;
    return dxdt;
}
// ----------------------------------------------------------
// analytic solution 
double sol(double x, double t){
    double solution;
    solution = exp(t);
    solution = 2.0 * exp(-0.5 * 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;
  int icount=0;
  // open files
  fp1 = fopen(filename1, "w");
  fp2 = fopen(filename2, "w");
  // initial set up
  t = T0;
  x = X0;
    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(x,t),x-sol(x,t));
          fprintf(fp1,"%12.5f %12.5f\n", t,x);
          fprintf(fp2,"%12.5f %12.5f\n", t,sol(x,t));
      }
      icount += 1;
    // Forward Difference
    x += dt * rhs(x,t);
    // next t
    t += dt;
  } // end of t-loop
  // close files 
  fclose(fp1);
  fclose(fp2);
}