PROGRAM d9r2(input,output);
(* driver for routine ZBRAC *)
VAR
   succes : boolean;
   i : integer;
   x1,x2 : real;

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

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

(*$I ZBRAC.PAS *)

BEGIN
   writeln('bracketing values:':21,'function values:':23);
   writeln('x1':6,'x2':10,
      'bessj0(x1)':21,'bessj0(x2)':12);
   FOR i := 1 to 10 DO BEGIN
      x1 := i;
      x2 := x1+1.0;
      zbrac(x1,x2,succes);
      IF  (succes)  THEN BEGIN
         writeln(x1:7:2,x2:10:2,' ':7,
            fx(x1):12:6,fx(x2):12:6)
      END
   END
END.
PROGRAM d9r3(input,output);
(* driver for routine ZBRAK *)
CONST
   n=100;
   nbmax=20;
   x1=1.0;
   x2=50.0;
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
VAR
   i,nb : integer;
   xb1,xb2 : glnbmax;

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

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

(*$I ZBRAK.PAS *)

BEGIN
   nb := nbmax;
   zbrak(x1,x2,n,xb1,xb2,nb);
   writeln;
   writeln('brackets for roots of bessj0:');
   writeln('lower':22,'upper':10,
         'f(lower)':16,'f(upper)':10);
   FOR i := 1 to nb DO BEGIN
      writeln('  root ',i:2,' ':4,xb1[i]:10:4,xb2[i]:10:4,
         ' ':4,fx(xb1[i]):10:4,fx(xb2[i]):10:4)
   END
END.
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.
PROGRAM d9r11(input,output);
(* driver for routine ZROOTS *)
CONST
   m=4;
   twomp2=10;   (* twomp2=(2*m+2) *)
TYPE
   glcarray = ARRAY [1..twomp2] OF real;
   gl2array = ARRAY [1..2] OF real;
VAR
   i : integer;
   polish : boolean;
   a,roots : glcarray;

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

(*$I ZROOTS.PAS *)

BEGIN
   a[1] := 0.0; a[2] := 2.0;
   a[3] := 0.0; a[4] := 0.0;
   a[5] := -1.0; a[6] := -2.0;
   a[7] := 0.0; a[8] := 0.0;
   a[9] := 1.0; a[10] := 0.0;
   writeln('Roots of the polynomial x^4-(1+2i)*x^2+2i');
   polish := false;
   zroots(a,m,roots,polish);
   writeln;
   writeln('Unpolished roots:');
   writeln('root #':14,'real':13,'imag.':13);
   FOR i := 1 to m DO BEGIN
      writeln(i:11,' ':5,roots[2*i-1]:12:6,roots[2*i]:12:6)
   END;
   writeln;
   writeln('Corrupted roots:');
   FOR i := 1 to m DO BEGIN
      roots[2*i-1] := roots[2*i-1]*(1+0.01*i);
      roots[2*i] := roots[2*i]*(1+0.01*i)
   END;
   writeln('root #':14,'real':13,'imag.':13);
   FOR i := 1 to m DO BEGIN
      writeln(i:11,' ':5,roots[2*i-1]:12:6,roots[2*i]:12:6)
   END;
   polish := true;
   zroots(a,m,roots,polish);
   writeln;
   writeln('Polished roots:');
   writeln('root #':14,'real':13,'imag.':13);
   FOR i := 1 to m DO BEGIN
      writeln(i:11,' ':5,roots[2*i-1]:12:6,roots[2*i]:12:6)
   END
END.
