#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 rk4(double y[],double dydx[],int &n,double &x,double &h,double yout [])
{int i;
double yt[11],dyt[11],dym[11],hh,h6,xh;
hh=h*0.5;
h6=h/6.0;
xh=x+hh;
for(i=1;i<=n;i++)
{yt[i]=y[i]+hh*dydx[i];
}
derivs(xh,yt,dyt);
for(i=1;i<=n;i++)
{yt[i]=y[i]+hh*dyt[i];
}
derivs(xh,yt,dym);
for(i=1;i<=n;i++)
{yt[i]=y[i]+h*dym[i];
dym[i]=dyt[i]+dym[i];
}
derivs(x+h,yt,dyt);
for(i=1;i<=n;i++)
{yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
}
for(i=1;i<=10;i++)
{dym[i]=0.0;
dyt[i]=0.0;
yt[i]=0.0;}
}


double bessj0(double x)
{double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
double y,bbb,ccc,aaa,temp,eee,ddd,ax,xx,z;
p1=1.0;                 p2=-0.001098628627;
p3=0.00002734510407;    p4=-0.00002073370639;
p5=2.093887211e-07;
q1=-0.01562499995;      q2=0.0001430488765;
q3=-0.000006911147651;  q4=7.621095161e-07;
q5=-9.34945152e-08;
r1=57568490574.0; r2=-13362590354.0;
r3=651619640.7; r4=-11214424.18;
r5=77392.33017; r6=-184.9052456;
s1=57568490411.0; s2=1029532985.0;
s3=9494680.718; s4=59272.64853;
s5=267.8532712; s6=1.0;
if(fabs(x)<8.0)
{y=x*x;
bbb=y*(r4+y*(r5+y*r6));
aaa=r1+y*(r2+y*(r3+bbb));
ccc=y*(s3+y*(s4+y*(s5+y*s6)));
temp=aaa/(s1+y*(s2+ccc));
}
else
{
	ax=fabs(x);
	z=8.0/ax;
	y=z*z;
	xx=ax-0.785398164;
	ccc=y*(p3+y*(p4+y*p5));
	aaa=p1+y*(p2+ccc);
	ddd=y*(q3+y*(q4+y*q5));
	eee=z*sin(xx)*(q1+y*(q2+ddd));
	temp=sqrt(0.636619772/ax)*(cos(xx)*aaa-eee);
}
return temp;
}


int sgn(double x)
{int temp;
if(x>0) temp=1;
if(x=0) temp=0;
if(x<0) temp=-1;
return temp;
}
double bessj1(double x)
{double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
double bbb,ccc,aaa,temp,ax,xx,z,y;
r1=72362614232.0; r2=-7895059235.0;
r3=242396853.1; r4=-2972611.439;
r5=15704.4826; r6=-30.16036606;
s1=144725228442.0;s2=2300535178.0;
s3=18583304.74; s4=99447.43394;
s5=376.9991397; s6=1.0;
p1=1.0; p2=0.00183105;
p3=-0.00003516396496; p4=0.000002457520174;
p5=-0.000000240337019; 
q1=0.04687499995; q2=-0.0002002690873;
q3=0.000008449199096;  q4=-0.00000088228987;
q5=0.000000105787412;
if(fabs(x)<8.0)
{y=x*x;
aaa=r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))));
bbb=s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))));
temp=x*aaa/bbb;
}
else
{ax=fabs(x);
z=8.0/ax;
y=z*z;
xx=ax-2.356194491;
aaa=p1+y*(p2+y*(p3+y*(p4+y*p5)));
bbb=q1+y*(q2+y*(q3+y*(q4+y*q5)));
ccc=sqrt(0.636619772/ax);
temp=ccc*(cos(xx)*aaa-z*sin(xx)*bbb*sgn(x));
}
return temp;
}



double bessj(int n, double x)
{int iacc,j,m;
double temp,bigno,bigni,ax,tox,bj,bjp,bjm,besj,jsum,sum;
iacc=40;
bigno=10000000000.0;
bigni=0.0000000001;
if(n<2)
{cout<<"bad argument n in bassj";
_c_exit();
}
ax=fabs(x);
if(ax==0)
{
	temp=0.0;
}
else if(ax>float(n))
{tox=2.0/ax;
bjm=bessj0(ax);
bj=bessj1(ax);
for(j=1;j<=n-1;j++)
{bjp=j*tox*bj-bjm;
bjm=bj;
bj=bjp;
}
temp=bj;
}
else
{tox=2.0/ax;
m=2*int(((n+int(sqrt(iacc*n))))/2);
besj=0.0;
jsum=0;
sum=0.0;
bjp=0.0;
bj=1.0;
for(j=m;j>=1;j--)
{bjm=j*tox*bj*-bjp;
bjp=bj;
bj=bjm;
if(fabs(bj)>bigno)
{bj=bj*bigni;
bjp=bjp*bigni;
besj=besj*bigni;
sum=sum*bigni;
}
if(jsum!=0) sum=sum+bj;
jsum=1-jsum;
if(j==n) besj=bjp;
}
sum=2.0*sum-bj;
temp=besj/sum;
}
return temp;
}

 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;
	}
 }