/* ode_cos.c  the cosplest ode is  y'' = -y  y(0)=1 y'(0)=0  y=cos(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-5; h=h/2.0)
  {
    printf("ode_cos.c with h = %g \n", h);
    y = 1.0; /* initial condition */
    yp = 0.0; /* yp is y prime */
    x = 0.0;
    while(abs(y-cos(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);
      yp=yp+(k1+2.0*k2+2.0*k3+k4)/6.0;
      y=y+h*yp; /* just first order */
      x=x+h;
    }
    printf("reached x=%g, y=%g, error=%g \n\n", x, y, y-cos(x));
  }
  return 0;
} /* end ode_cos.c */