/* Runge Kutta algorithm for first-order differential equations*/

#include <stdio.h>

#define dist 0.5                        /* stepsize in t */ 
#define MAX 4                           /* max for t */

FILE *output;                           /* internal filename */
 
main()
{
double t, y;
int j;
 
double runge4(double x, double y, double step); /* Runge-Kutta function */

double f(double x, double y);           /* function for derivative */

output=fopen ("runge.dat", "w");        /* external filename */

y=1;                                    /* initial position */
fprintf(output, "0\t%f\n", y);

for (j=1; j*dist<=MAX; j++)             /* the time loop */
{
   t=j*dist;
   y=runge4(t, y, dist);

   fprintf(output, "%f\t%f\n", t, y);
}

fclose(output);
}

double runge4(double x, double y, double step)
{
double h=step/2.0, k1, k2, k3, k4;

k1=step*f(x, y);
k2=step*f(x+h, y+k1/2.0);
k3=step*f(x+h, y+k2/2.0);
k4=step*f(x+step, y+k3);

return(y+(k1+2*k2+2*k3+k4)/6.0);
}

double  f(double x, double y)
{
return(-y);
}