PROGRAM d9r7(input,output);
(* driver for routine ZBRENT *)
CONST
   n=100;
   nbmax=20;
   x1=1.0;
   x2=50.0;
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
VAR
   i,nb : integer;
   root,tol : real;
   xb1,xb2 : glnbmax;

(*$I MODFILE.PAS *)
(*$I BESSJ0.PAS *)

FUNCTION fx(x: real): real;
BEGIN
   fx := bessj0(x)
END;

(*$I ZBRAK.PAS *)

(*$I ZBRENT.PAS *)

BEGIN
   nb := nbmax;
   zbrak(x1,x2,n,xb1,xb2,nb);
   writeln;
   writeln('Roots of bessj0:');
   writeln('x':18,'f(x)':15);
   FOR i := 1 to nb DO BEGIN
      tol := (1.0e-6)*(xb1[i]+xb2[i])/2.0;
      root := zbrent(xb1[i],xb2[i],tol);
      writeln('root ',i:2,'  ',root:12:6,fx(root):14:6)
   END
END.