/* 
************************************************************************
*  runge.c: 4th order Runge-Kutta solution for an ordinary differential
*  	  equation
************************************************************************
*/
       
#include <stdio.h>
#include <math.h>

#define MIN 0.0                               /* minimum x */
#define MAX 5.0                              /* maximum x */
#define PI 3.1415926
 
main()
{
   double f(double x, double y);
   
   double x, y, dx, h, dy, yreal, yerror;

   double t1, t2, t3,                /* temporary storage */
          k1, k2, k3, k4;          /* for Runge-Kutta  */
   
   FILE *output;                              /* save data in rk4.dat */
   output = fopen("runge.dat","w");

   printf("Input: dx\n");
   scanf("%lf", &dx);
   printf("dx=%g\n", dx);
   
   x = MIN;                                /* initial position  */
   y = 0.0;
   yreal = y;
   yerror = 0.0;
   
   fprintf(output, "%f\t%f\t%f\t%e\n", x, y, yreal, yerror);

   while(x <= MAX)
   {
      h=dx/2.0;
      k1 = dx*f(x, y);
      t1 = y+0.5*k1;
      k2 = dx*f(x+h, t1);
      t2 = y+0.5*k2;
      k3 = dx*f(x+h, t2);
      t3 = y+k3;
      k4 = dx*f(x + dx, t3);
      x += dx;
      y += (k1+2*k2+2*k3+k4)/6.0;      

      yreal=sin(2.0*PI*x);
      yerror=y-yreal;      
      fprintf(output, "%f\t%f\t%f\t%e\n", x, y, yreal, yerror);
   }
   printf("data stored in runge.dat\n");
   fclose(output);
}
/*-----------------------end of main program--------------------------*/

/* definition of equations  */
double  f(double x, double y)
{
   return(2.0*PI*cos(2.0*PI*x));               /* RHS of first equation */
}