#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#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 d14r3
    //driver for routine rkqc
	int n,i;
	double y[5], dydx[5], yscal[5],x,htry,eps,hdid,hnext;
    n = 4;  
    x = 1.0;
	hnext=0.0;
	hdid=0.0;
	eps=0.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));
        rkqc(y, dydx, n, x, htry, eps, yscal, hdid, hnext);
        cout<<setprecision(7)<<setw(12)<<eps;
        cout<<setprecision(6)<<setw(13)<<htry;
        cout<<setprecision(6)<<setw(13)<<hdid;
        cout<<setprecision(6)<<setw(13)<<hnext<<endl;
	}
}