PROGRAM d15r6(input,output);
(* driver for routine BSSTEP *)
CONST
   n=4;
   glimax=11;
   glnmax=10;
   glncol=7;
TYPE
   glyarray = ARRAY [1..n] OF real;
   glnarray = glyarray;
VAR
   glx : ARRAY [1..glimax] OF real;
   gld : ARRAY [1..glnmax,1..glncol] OF real;
   eps,hdid,hnext,htry,x : real;
   i : integer;
   y,dydx,yscal : glyarray;

PROCEDURE derivs(x: real; y: glyarray; VAR dydx: glyarray);
(* Programs using DERIVS must define the type
TYPE
   glyarray = ARRAY [1..4] OF real;
in the main routine. *)
BEGIN
   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]
END;

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

(*$I BESSJ1.PAS *)

(*$I BESSJ.PAS *)

(*$I MMID.PAS *)

(*$I RZEXTR.PAS *)

(*$I BSSTEP.PAS *)

BEGIN
   x := 1.0;
   y[1] := bessj0(x);
   y[2] := bessj1(x);
   y[3] := bessj(2,x);
   y[4] := bessj(3,x);
   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];
   FOR i := 1 to n DO BEGIN
      yscal[i] := 1.0
   END;
   htry := 1.0;
   writeln;
   writeln('eps':10,'htry':12,'hdid':12,'hnext':12);
   FOR i := 1 to 15 DO BEGIN
      eps := exp(-i);
      bsstep(y,dydx,n,x,htry,eps,yscal,hdid,hnext);
      writeln('   ',eps:11,htry:8:2,hdid:14:6,hnext:12:6)
   END
END.