OK. So here is the main code.

Thank you.

Code:

#include <iostream>
#include <iomanip>
#include "nr.h"
using namespace std;
using std::fixed;
// Driver for routine odeint
DP dxsav; // defining declarations
int kmax,kount;
Vec_DP *xp_p;
Mat_DP *yp_p;
int nrhs; // counts function evaluations
void derivs(const DP x, Vec_I_DP &y, Vec_O_DP &dydx)
{
nrhs++;
dydx[0] = -y[1];
dydx[1] = y[0] - (1.0/x)*y[1];
dydx[2] = y[1] - (2.0/x)*y[2];
dydx[3] = y[2] - (3.0/x)*y[3];
}
int main(void)
{
const int N = 4, KMAX = 100;
int i, nbad, nok;
DP eps = 1.0e-4, h1 = 0.1, hmin = 0.0, x1 = 1.0, x2 = 10.0;
Vec_DP ystart(N);
ystart[0] = NR::bessj0(x1);
ystart[1] = NR::bessj1(x1);
ystart[2] = NR::bessj(2,x1);
ystart[3] = NR::bessj(3,x1);
nrhs = 0;
dxsav = (x2-x1)/20.0;
kmax = KMAX;
xp_p = new Vec_DP(KMAX);
yp_p = new Mat_DP(N,KMAX);
Vec_DP &xp = *xp_p;
Mat_DP &yp = *yp_p;
NR::odeint(ystart,x1,x2,eps,h1,hmin,nok,nbad,derivs,NR::rkqs);
cout << endl << "successful steps:" << setw(14) << " ";
cout << setw(4) << nok << endl;
cout << "bad steps:" << setw(21) << " " << setw(4) << nbad << endl;
cout << "function evaluations:" << setw(10) << " ";
cout << setw(4) << nrhs << endl;
cout << endl << "stored intermediate values: ";
cout << setw(4) << kount << endl;
cout << endl << setw(8) << "x" << setw(18) << "integral";
cout << setw(16) << "bessj(3,x)" << endl;
// cout << fixed << setprecision(6);
for (i=0;i<kount;i++) {
cout << setw(10) << xp[i] << setw(16) << yp[3][i];
cout << setw(15) << NR::bessj(3,xp[i]) << endl;
}
delete yp_p;
delete xp_p;
return 0;
}