void bsstep(double y[],double dydx[],int nv,double& x1,double htry,double eps,double yscal[],double& hdid,double &hnext)
{
	int imax,nuse,i,j;
	double h,xsav,errmax;
	double one,shrink,grow,xest;
	double yerr[11],ysav[11],dysav[11],yseq[11],nseq[12];
    imax = 11;
    nuse = 7;
    one = 1.0;
    shrink = 0.95;
    grow = 1.2;   
    nseq[1] = 2.0;   nseq[2] = 4.0;  nseq[3] = 6.0;  nseq[4] = 8.0;
	nseq[5] = 12.0;  nseq[6] = 16.0; nseq[7] = 24.0; nseq[8] = 32.0;
	nseq[9] = 48.0;  nseq[10]= 64.0; nseq[11]= 96.0;
    h = htry;
    xsav = x1;
    for (i = 1; i<=nv; i++)
	{
        ysav[i] = y[i];
        dysav[i] = dydx[i];
    }
    do
	{
        for (i = 1; i<=imax; i++)
		{
            mmid(ysav,dysav,nv,xsav,h,nseq[i],yseq);
		
            xest = (h / nseq[i])* (h / nseq[i]);
            rzextr(i,xest,yseq,y,yerr,nv,nuse);
            if (i > 3)
			{
                errmax = 0.0;
                for (j = 1; j<=nv; j++)
				{
                    if (fabs(yerr[j] / yscal[j]) > errmax)
					{
                        errmax = fabs(yerr[j] / yscal[j]);
                    }
                }
                errmax = errmax / eps;
                if (errmax < one)
				{
                    x1 = x1 + h;
                    hdid = h;
                    if (i == nuse)
                        hnext = h * shrink;
                    else if (i == (nuse - 1))
                        hnext = h * grow;
                    else
                        hnext = (h * nseq[nuse - 1]) / nseq[i];
                    for (i=1; i<=12; i++)
					{
                     nseq[i]=0.0;
					}
					 for (i=1; i<=11; i++)
					 {
					    yseq[i]=0.0;
						dysav[i]=0.0;
						ysav[i]=0.0;
						yerr[i]=0.0;
					}
                   return;
                }
			}
        }
        h = 0.25 * h / pow(2 , ((imax - nuse) / 2));
	}
    while ((x1 + h)!= x1);
    cout<< " step size underflow"<<endl;
}