#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 d14r5
    //driver for routine mmid
	int i;
	double y[5], yout[5], dydx[5],nvar,x1,htot,b1,b2,b3,b4,xf;
    nvar = 4;
    x1 = 1.0;
    htot = 0.5;
    y[1] = bessj0(x1);
    y[2] = bessj1(x1);
    y[3] = bessj(2, x1);
    y[4] = bessj(3, x1);
    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];
    xf = x1 + htot;
    b1 = bessj0(xf);
    b2 = bessj1(xf);
    b3 = bessj(2, xf);
    b4 = bessj(3, xf);
    cout<<"First four Bessl functions"<<endl;
    for( i = 5; i<=50; i+=5)
	{
    mmid(y, dydx, nvar, x1, htot, i, yout);
    cout<<"x= "<<x1<<" to "<<x1 + htot<<" in "<<i<<" steps"<<endl;
    cout<<"    Integration      Bessj"<<endl;
    cout<<setw(13)<<yout[1];    cout<<setw(15)<<b1<<endl;
    cout<<setw(13)<<yout[2];    cout<<setw(15)<<b2<<endl;
    cout<<setw(13)<<yout[3];    cout<<setw(15)<<b3<<endl;
    cout<<setw(13)<<yout[4];    cout<<setw(15)<<b4<<endl;
    }
}