/* vsrk4.c: Fourth-order Runge-kutta method with variable step. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
extern double f( double x, double y );
double rk4( double x, double y, double hh );
int vsrk4(double hm[2],double x01[2],double x[],double y[],double eps)
{
int i;
double xx,h,hmax,hmin,x0,x1,yn,ynh,delta;
hmin=hm[0]; hmax=hm[1];
x0=x01[0];  x1=x01[1];
h=hmax;
xx=x0;
i=0;
x[0]=xx;
while( xx < x1 )
  {
again:  if( xx+h > x1 ) h=x1-xx;
  yn=rk4(xx,y[i],h);
  ynh=rk4(xx,y[i],h/2.0);
  ynh=rk4(xx+h/2.0,ynh,h/2.0);
  delta=fabs((yn-ynh)/15.0);
  if( ( delta < eps )||( h == hmin ) )
    {
    xx=xx+h;
    i=i+1;
    x[i]=xx;
    y[i]=ynh;
    if( i >= 99 ) goto endd;
    if( delta < 0.05*eps )
      {
      h=2.0*h;
      if( h > hmax ) h=hmax;
      }
    }
  else
    {
    h=h/2.0;
    if( h < hmin ) h=hmin;
    goto again;
    }
  }
endd:return(i);
}
/**************** fourth order Runge-Kutta *************************/
double rk4( double xn, double yn, double h )
{
double ynp,k1,k2,k3,k4;
k1=h*f(xn,yn);
k2=h*f(xn+0.5*h,yn+0.5*k1);
k3=h*f(xn+0.5*h,yn+0.5*k2);
k4=h*f(xn+h,yn+k3);
ynp=yn+(k1+2.0*k2+2.0*k3+k4)/6.0;
return(ynp);
}