void rk4(y,dydx,x,h,yout,derivs)
double y[],dydx[],x,h,yout[];
void (*derivs)();	/* ANSI: void (*derivs)(double,double *,double *); */
{
	double hh,h6,dym[6],dyt[6],yt[10];
	int xh;
	extern double Cinc();

	hh=h*0.5;
	h6=h/6.0;
	xh=x+Cinc(hh);
	yt[6]=y[6]; yt[7]=y[7]; yt[8]=y[8]; yt[9]=y[9];

	yt[0]=y[0]+hh*dydx[0];
	yt[1]=y[1]+hh*dydx[1];
	yt[2]=y[2]+hh*dydx[2];
	yt[3]=y[3]+hh*dydx[3];
	yt[4]=y[4]+hh*dydx[4];
	yt[5]=y[5]+hh*dydx[5];
	(*derivs)(dyt, xh,yt);

	yt[0]=y[0]+hh*dyt[0];
	yt[1]=y[1]+hh*dyt[1];
	yt[2]=y[2]+hh*dyt[2];
	yt[3]=y[3]+hh*dyt[3];
	yt[4]=y[4]+hh*dyt[4];
	yt[5]=y[5]+hh*dyt[5];
	(*derivs)(dym, xh,yt);

	yt[0]=y[0]+h*dym[0];  dym[0] += dyt[0];
	yt[1]=y[1]+h*dym[1];  dym[1] += dyt[1];
	yt[2]=y[2]+h*dym[2];  dym[2] += dyt[2];
	yt[3]=y[3]+h*dym[3];  dym[3] += dyt[3];
	yt[4]=y[4]+h*dym[4];  dym[4] += dyt[4];
	yt[5]=y[5]+h*dym[5];  dym[5] += dyt[5];
	(*derivs)(dyt, (int)(x+Cinc(h)),yt);

	yout[0]=y[0]+h6*(dydx[0]+dyt[0]+2.0*dym[0]);
	yout[1]=y[1]+h6*(dydx[1]+dyt[1]+2.0*dym[1]);
	yout[2]=y[2]+h6*(dydx[2]+dyt[2]+2.0*dym[2]);
	yout[3]=y[3]+h6*(dydx[3]+dyt[3]+2.0*dym[3]);
	yout[4]=y[4]+h6*(dydx[4]+dyt[4]+2.0*dym[4]);
	yout[5]=y[5]+h6*(dydx[5]+dyt[5]+2.0*dym[5]);
}