// Solve 2nd order differential equation using Euler method
//                                2019/1   Hisaaki Shinkai
// compile as ::   gcc -lm -o DE2.exe DE2.c
// execute as ::   ./DE2.exe
// output files :: output.analytic     x  y
//                 output.numerical    x  y
// plot graphs using gnuplot
//      plot "output.analytic" with lines, "output.numerical" with lines
#include <stdio.h>
#include <math.h>

#define Y0 1.0       /* Initial Value  y(0) */
#define DY0 0.0      /* Initial Value  y'(0) */
#define T0 0.0       /* starting time t0 */
#define TMAX 10.0    /* ending time tmax */

int main(void)
{
  char filename1[] = "output.numerical";
  char filename2[] = "output.analytic";
  FILE *fp1, *fp2;
  double t=0.0, dt=0.05;
  double y=0.0;
  double z=0.0;
  double dydt=0.0, d2ydt2=0.0;
  // Opening message
  printf("Welcome to DE2.c \n");
  // open files
  fp1 = fopen(filename1, "w");
  fp2 = fopen(filename2, "w");
  if (fp1 == NULL)
    {
      printf("opening file error %s\n", filename1);
	  return(-1);
    }
	if (fp2 == NULL)
    {
      printf("opening file error %s\n", filename2);
	  return(-1);
    }
  // t-Loop
  t = T0;
  y = Y0;
  dydt = DY0;
  while(t < TMAX){
    // Set your problem below
    //     d2ydt2 = right-hand side of the 2nd order DE
    d2ydt2 = -3.0 * dydt -2.0 *y + 3.0 * sin(2.0 * t);
    d2ydt2 = -0.1 * y + 0.3 * sin(t);
    d2ydt2 = -0.1 * y + 0.3 * sin(t) * sin(t);
    d2ydt2 = -4.0 * dydt - 4.0 * y;
    d2ydt2 = -2.0 * dydt - 2.0 * y + 3.0;
    d2ydt2 = -4.0 * y;
    //     Set your problem  ... end
    // Write your analytic solution below
    z = y;
    z = exp(-0.2 * t);
    z = 1.0 * cos(2.0 * t);
    //     Analytic solution ... end
    // output
	printf("%10.3f %11.5f %11.5f %12.8f \n", t,y,z,y-z);
	fprintf(fp1,"%12.5f %12.5f\n", t,y);
	fprintf(fp2,"%12.5f %12.5f\n", t,z);
    // Forward Difference
    dydt += d2ydt2 * dt;
    y += dydt * dt;
	// next t
    t += dt;
  } // end of t-loop
  // close files 
  fclose(fp1);
  fclose(fp2);
  return (0);
}