/* ode_exp.c  the simplest ode is  y' = y  y(0)=1  y=exp(x)  */
/*            my basic ode solver is Runge Kutta 4th         */
/*            experiment with h to see how far error < 0.001 */
#include <stdio.h>
#include <math.h>
#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))

/* for this ode we still need a function for y' in terms of x and y */
double f(double x, double y)
{
  return y; /* x is not needed, yet we use the standard function */
}

int main(int argc, char *argv[])
{
  double h, y, x, yp;
  double k1,k2,k3,k4;

  for(h=0.25; h>1.0e-6; h=h/2.0)
  {
    printf("ode_exp.c with h = %g \n", h);
    y = 1.0; /* initial condition */
    x = 0.0;
    while(abs(y-exp(x))<0.001) /* loop until error > 0.001 */
    {
      k1=h*f(x,y);
      k2=h*f(x+h/2.0,y+k1/2.0);
      k3=h*f(x+h/2.0,y+k2/2.0);
      k4=h*f(x+h,y+k3);
      y=y+(k1+2.0*k2+2.0*k3+k4)/6.0;
      x=x+h;
    }
    printf("reached x=%g, y=%g, error=%g \n\n", x, y, y-exp(x));
  }
  return 0;
} /* end ode_exp.c */