#include<math.h>
#include<iomanip.h>
#include<iostream.h>
#include<process.h>

int sgn(double& w)
{
	int t;
	if (w>0) t=1;
    if (w==0) t=0;
	if (w>0) t=-1;
	return t;

}
double bessj0(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double y,bbb,ccc,aaa,temp,eee,ddd,ax,xx,z;
    p1 = 1.0;                  p2 = -0.001098628627;
    p3 = 0.00002734510407;    p4 = -0.000002073370639;
    p5 = 2.093887211e-07;
    q1 = -0.01562499995;      q2 = 0.0001430488765;
    q3 = -0.000006911147651;  q4 = 7.621095161e-07;
    q5 = -9.34945152e-08;
    r1 = 57568490574.0;        r2 = -13362590354.0;
    r3 = 651619640.7;         r4 = -11214424.18;
    r5 = 77392.33017;         r6 = -184.9052456;
    s1 = 57568490411.0;        s2 = 1029532985.0;
    s3 = 9494680.718;         s4 = 59272.64853;
    s5 = 267.8532712;        s6 = 1.0;
    if (fabs(x) < 8.0)
	{
       y = x * x;
       bbb = y * (r4 + y * (r5 + y * r6));
       aaa = r1 + y * (r2 + y * (r3 + bbb));
       ccc = y * (s3 + y * (s4 + y * (s5 + y * s6)));
       temp = aaa / (s1 + y * (s2 + ccc));
	}
    else
	{
       ax = fabs(x);
       z = 8.0 / ax;
       y = z * z;
       xx = ax - 0.785398164;
       ccc = y * (p3 + y * (p4 + y * p5));
       aaa = p1 + y * (p2 + ccc);
       ddd = y * (q3 + y * (q4 + y * q5));
       eee = z * sin(xx) * (q1 + y * (q2 + ddd));
       temp = sqrt(0.636619772 / ax) * (cos(xx) * aaa - eee);
	}
	return temp;
}

  double bessj1(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double bbb,ccc,aaa,temp,ax,xx,z,y;
    r1 = 72362614232.0;           r2 = -7895059235.0;
    r3 = 242396853.1;            r4 = -2972611.439;
    r5 = 15704.4826;             r6 = -30.16036606;
    s1 = 144725228442.0;          s2 = 2300535178.0;
    s3 = 18583304.74;            s4 = 99447.43394;
    s5 = 376.9991397;            s6 = 1.0;
    p1 = 1.0;                     p2 = 0.00183105;
    p3 = -0.00003516396496;      p4 = 0.000002457520174;
    p5 = -0.000000240337019;
    q1 = 0.04687499995;          q2 = -0.0002002690873;
    q3 = 0.000008449199096;      q4 = -0.00000088228987;
    q5 = 0.000000105787412;
    if (fabs(x) < 8.0 )
	{
        y = x*x;
        aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
        bbb = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
        temp = x * aaa / bbb;
	}
    else
	{
        ax = fabs(x);
        z = 8.0 / ax;
        y = z*z;;
        xx = ax - 2.356194491;
        aaa = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
        bbb = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
        ccc = sqrt(0.636619772 / ax);
        temp = ccc * (cos(xx) * aaa - z * sin(xx) * bbb * sgn(x));
	}  
	return temp;
}
double bessj(int n, double x)
{
	int iacc,j,m;
	double temp,bigno,bigni,ax,tox,bj,bjp,bjm,besj,jsum,sum;
    iacc = 40;
    bigno = 10000000000.0;
    bigni = 0.0000000001;
    if (n < 2 )
	{
        cout<<"bad argument n in bassj";
        _c_exit();
    }
    ax = fabs(x);
    if( ax == 0) 
	{
        temp = 0.0;
	}
    else if (ax > float(n))
	{
        tox = 2.0 / ax;
        bjm = bessj0(ax);
        bj = bessj1(ax);
        for( j = 1;j<=n-1;j++)
		{
            bjp = j * tox * bj - bjm;
            bjm = bj;
            bj = bjp;
        }
        temp = bj;
	}
    else
	{
        tox = 2.0 / ax;
        m = 2 * int(((n + int(sqrt(iacc * n)))) / 2);
        besj = 0.0;
        jsum = 0;
        sum = 0.0;
        bjp = 0.0;
        bj = 1.0;
        for (j = m;j>=1;j--)
		{
            bjm = j * tox * bj - bjp;
            bjp = bj;
            bj = bjm;
            if (fabs(bj) > bigno)
			{
                bj = bj * bigni;
                bjp = bjp * bigni;
                besj = besj * bigni;
                sum = sum * bigni;
            }
            if (jsum != 0) sum = sum + bj;
            jsum = 1 - jsum;
            if (j == n) besj = bjp;
        }
        sum = 2.0 * sum - bj;
        temp = besj / sum;		
    }
	return temp;
}



void rzextr(int iest,double xest,double  yest[],double yz[],double  dy[],int nv,int nuse)
{
	int j,m1,k;
    double  fx[8],yy,v,c,b1,b,ddy;
    static double x[12],d[11][ 8];
    x[iest] = xest;
    if (iest == 1)
	{
        for (j = 1;j<= nv;j++)
		{
            yz[j] = yest[j];
            d[j][1] = yest[j];
            dy[j] = yest[j];
		}
	}
    else
	{
        m1 = iest;
        if (nuse < iest) m1 = nuse;
        for( k = 1; k<= m1 - 1;k++)
		{
            fx[k + 1] = x[iest - k] / xest;
        }
        for (j = 1;j<= nv;j++)
		{
            yy = yest[j];
            v = d[j][1];
            c = yy;
            d[j][1] = yy;
            for (k = 2;k<=m1;k++)
			{
                b1 = fx[k] * v;
                b = b1 - c;
                if( b != 0.0)
				{
                    b = (c - v) / b;
                    ddy = c * b;
                    c = b1 * b;
				}
                else
				{
                    ddy = v;
                }
                if (k != m1)  v = d[j][k];
                d[j][k] = ddy;
                yy = yy + ddy;
            }
            dy[j] = ddy;
            yz[j] = yy;
		}
    }
	for (j = 1;j<=7;j++)
	{
    fx[j]=0.0;
	}
    
}

void derivs(double x,double y[],double dydx[])
{
    dydx[1] = -y[2];
    dydx[2] = y[1] - (1.0 / x) * y[2];
    dydx[3] = y[2] - (2.0 / x) * y[3];
    dydx[4] = y[3] - (3.0 / x) * y[4];
}

void mmid(double y[],double dydx[], int &nvar, double &xs, double &htot, double 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]);
    }
}

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;
}

void main()
{
    	int n,i;
    n = 4;
    double  y[5],dydx[5],yscal[5],x;
	double htry,eps,hdid,hnext;
    x = 1.0;
    y[1] = bessj0(x);
    y[2] = bessj1(x);
    y[3] = bessj(2,x);
    y[4] = bessj(3,x);
    dydx[1] = -y[2];
    dydx[2] = y[1] - y[2];
    dydx[3] = y[2] - 2.0 * y[3];
    dydx[4] = y[3] - 3.0 * y[4];
    for (i = 1; i<=n; i++)
	{
        yscal[i] = 1.0;
    }
    htry = 0.1;
    cout<<endl;
    cout<< "      eps       htry        hdid          hnext"<<endl;
	cout<<setiosflags(ios::fixed);
    for (i = 1; i<=15; i++)
	{
        eps = exp(-float(i));
        bsstep(y,dydx,n,x,htry,eps,yscal,hdid,hnext);
		
        cout<<setprecision(7)<<setw(12)<<eps;
        cout<<setprecision(2)<<setw(8)<<htry;
        cout<<setprecision(6)<<setw(14)<<hdid;
        cout<<setprecision(6)<<setw(14)<<hnext<<endl;
    }
}