// Solve 1st order differential equation using Euler method
//                                2019/1   Hisaaki Shinkai
// compile as ::   gcc -lm -o DE1.exe DE1.c
// execute as ::   ./DE1.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 2.0       /* Initial Value  y(0) */
#define X0 0.0       /* starting coord x0 */
#define XMAX 10.0    /* ending coord xmax */

int main(void)
{
  char filename1[] = "output.numerical";
  char filename2[] = "output.analytic";
  FILE *fp1, *fp2;
  double x=0.0, dx=0.01;
  double y=0.0;
  double z=0.0;
  double dydx=0.0;
  // Opening message
  printf("Welcome to DE1.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);
    }
  // x-Loop
  x = X0;
  y = Y0;
  while(x < XMAX){
    // *** Set your problem below
    //     dydx = right-hand side of the 1st order DE
    dydx = -0.5 * cos(x) * y;
    dydx = -0.5 * y + exp(-0.2 * x);
    dydx = -0.5 * y + sin(x);
    dydx = -0.5 * y + 1.0;
    dydx = -0.2 * x * y;
    dydx = -0.2 * y;
    //     Set your problem  ... end
    // *** Write your analytic solution below
    z = exp(x);
    z = 2.0 * exp(-0.2 * x);
    //     Analytic solution ... end
    // output
	printf("%10.3f %11.5f %11.5f %12.8f \n", x,y,z,y-z);
	fprintf(fp1,"%12.5f %12.5f\n", x,y);
	fprintf(fp2,"%12.5f %12.5f\n", x,z);
    // Forward Difference
    y += dydx * dx;
    // next x
    x += dx;
  } // end of x-loop
  // close files
  fclose(fp1);
  fclose(fp2);
  return (0);
}