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

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 main()
{
    //program d14r6
    //driver for routine bsstep
	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;
    }
}