void rk4(double y[], double dydx[], int &n, double &x, double &h, double yout[])
{
	int i;
    double yt[11], dyt[11], dym[11],hh,h6,xh;
    hh = h * 0.5;
    h6 = h / 6.0;
    xh = x + hh;
    for( i = 1; i<=n; i++)
	{
        yt[i] = y[i] + hh * dydx[i];
    }
    derivs(xh, yt, dyt);
    for( i = 1; i<=n; i++)
	{
        yt[i] = y[i] + hh * dyt[i];
    }
    derivs(xh, yt, dym);
    for( i = 1; i<=n; i++)
	{
        yt[i] = y[i] + h * dym[i];
        dym[i] = dyt[i] + dym[i];
    }
    derivs(x + h, yt, dyt);
    for (i = 1; i<=n; i++)
	{
        yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
    }
	for(i=1; i<=10; i++)
	{
		dym[i]=0.0;
		dyt[i]=0.0;
		yt[i]=0.0;
	}
}