#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 d14r1
    //driver for routine rk4
	int n,i,j;
	double y[5], dydx[5], yout[5],h,x;
    n = 4;
    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];
    cout<<"Bessel function:  j0      j1      j2      j3"<<endl;
    for (i = 1; i<=5; i++)
	{
        h = 0.2 * i;
        rk4(y, dydx, n, x, h, yout);
        cout<<endl;
        cout<<"for a step size of: "<<h<<endl;
        cout<<"  rk4:    ";
        for (j = 1; j<=4; j++)
		{
            cout<<setw(12)<<yout[j];
        }
		cout<<endl;
        cout<<"  actual:  ",
        cout<<setw(12)<<bessj0(x + h);
        cout<<setw(12)<<bessj1(x + h);
        cout<<setw(12)<<bessj(2, x + h);
        cout<<setw(12)<<bessj(3, x + h)<<endl;
    }
}