//Linda Cobb

//adapted from 'numerical recipes in c'



//http://www.timestocome.com





#define NRANSI

#include "nrutil.h"



void rk4(double y[], double  dydx[], int n, double x, double h, double yout[],

	void (*derivs)(double, double [], double []))

{



	int i;

	double xPlush, halfstep, h6, *dym, *dyt, *yt;



	dym = dvector(1,n);

	dyt = dvector(1,n);

	yt  = dvector(1,n);

	halfstep  = h*0.5;

	h6  = h/6.0;

	



	for (i=1;i<=n;i++)                       //for each ode do an euler step

	  yt[i] = y[i] + halfstep*dydx[i];       //dydx given for here



	(*derivs)(x + halfstep, yt, dyt);        //2nd rk4 step  

	for (i=1;i<=n;i++) 

	  yt[i] = y[i] + halfstep*dyt[i];



	(*derivs)(x + halfstep, yt, dym);        //3rd step

	for (i=1;i<=n;i++) {

	  yt[i] = y[i] + h*dym[i];

	  dym[i] += dyt[i];

	}



	(*derivs)(x + h, yt, dyt);                 //4th step

	for (i=1;i<=n;i++)                         //weight total

		yout[i] = y[i] + h6*(dydx[i] + dyt[i] + 2.0*dym[i]);





	free_dvector(yt,1,n);

	free_dvector(dyt,1,n);

	free_dvector(dym,1,n);

}

#undef NRANSI