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

double func(double x)
{
	double t;
	t=bessj0( x);
	return t;
}

void  funcd(double x, double &fn,double& df)
{
    fn = bessj0(x);
    df = -bessj1(x);
}

void main()
{
    //program d10r9
    //driver for routine rtsafe
	int n,nbmax,nb,i;
	double  x1,x2,xacc,root,xb1[21], xb2[21],ttt;
    n = 100;
    nbmax = 20;
    x1 = 1.0;
    x2 = 50.0;
    nb = nbmax;
    zbrak(x1, x2, n, xb1, xb2, nb);
    cout<<endl;
    cout<<"Roots of Bessj0:"<<endl;
    cout<<endl;
    cout<< "                    x              f(x)"<<endl;
    for (i = 1; i<=nb; i++)
	{
        xacc = (0.0000010) * (xb1[i] + xb2[i]) / 2.0;
        root =rtsafe(xb1[i], xb2[i], xacc);
        cout<< setw(8)<<"root "<<i;
        cout<< setw(14)<<root;
		ttt=bessj0(root);
        cout<< setw(19)<<ttt<<endl;
    }
}