#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 d14r4
    //driver for routine odeint
	int i,nvar;
	double vstart[5], xp[201], yp[10][201],x1,x2,eps,h1,hmin,nok,nbad;
    nvar = 4;
    x1 = 1.0;
    x2 = 10.0;
    vstart[1] = bessj0(x1);
    vstart[2] = bessj1(x1);
    vstart[3] = bessj(2, x1);
    vstart[4] = bessj(3, x1);
    eps = 0.0001;
    h1 = 0.1;
    hmin = 0.0;
    kmax = 100;
    dxsav = (x2 - x1) / 20.0;
    odeint(vstart, nvar, x1, x2, eps, h1, hmin, nok, nbad, xp, yp);
    cout<<endl;
    cout<<"Successful step: "<<nok<<endl;
    cout<<"Bad step: "<<nbad<<endl;
    cout<<"Stored intermediate values: "<<kount<<endl;
    cout<<endl;
    cout<<"       x         Integral      Bessj(3,x)"<<endl;
    for( i = 1; i<=kount; i++)
	{
        cout<<setprecision(4)<<setw(10)<<xp[i];
        cout<<setprecision(6)<<setw(15)<<yp[4][i];
        cout<<setprecision(6)<<setw(15)<<bessj(3, xp[i])<<endl;
    }
}