#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>

int cint(double x)
{
	int temp;
	double iprt;
	if (x>0)
	{
	x=modf(x,&iprt);
		if(fabs(x)<0.5)
			temp=int(iprt);
		else
			temp=int(iprt+1);
	}
	else if(x==0)
		temp=0;
	else
	{
		x=modf(x,&iprt);
		if(fabs(x)<0.5)
			temp=int(iprt);
		else
			temp=int(iprt)-1;
	}
		return temp;
}

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 d14r2
    //driver for routine rkdumb
	int i,j;
	double vstart[5], xx[201], y[11][201],nvar,nstep,x1,x2;
    nvar = 4;
    nstep = 150;
    x1 = 1.0;
    vstart[1] = bessj0(x1);
    vstart[2] = bessj1(x1);
    vstart[3] = bessj(2,x1);
    vstart[4] = bessj(3,x1);
    x2 = 20.0;
    rkdumb(vstart, nvar, x1, x2, nstep, xx, y);
    cout<<endl;
    cout<<"    x          Integrated        Bessj3"<<endl;
    for( i = 1; i<=cint(nstep / 10); i++)
	{
        j = 10 * i;
        cout<<setw(8)<<xx[j];
        cout<<setw(16)<<y[4][j];
        cout<<setw(16)<<bessj(3, xx[j])<<endl;
    }
}