void mmid(double y[],double dydx[], double &nvar, double &xs, double &htot, int nstep, double yout[])
{
	int n,i;
    double ym[11],yn[11],x,swap,h,h2;
    h = htot / nstep;
    for( i = 1; i<=nvar; i++)
	{
        ym[i] = y[i];
        yn[i] = y[i] + h * dydx[i];
    }
    x = xs + h;
    derivs(x, yn, yout);
    h2 = 2.0 * h;
    for( n = 2; n<=nstep; n++)
	{
        for (i = 1; i<=nvar; i++)
		{
            swap = ym[i] + h2 * yout[i];
            ym[i] = yn[i];
            yn[i] = swap;
        }
        x = x + h;
        derivs(x, yn, yout);
    }
    for( i = 1; i<=nvar; i++)
	{
        yout[i] = 0.5 * (ym[i] + yn[i] + h * yout[i]);
    }
}