PROGRAM d15r6(input,output);
(* driver for routine BSSTEP *)
   
(*$I MODFILE.PAS *)

CONST
   n = 4;
   RzextrImax = 11;
   RzextrNmax = 10;
   RzextrNcol = 7;
TYPE
   RealArrayNVAR = ARRAY [1..n] OF real;
VAR
   RzextrX: ARRAY [1..RzextrImax] OF real;
   RzextrD: ARRAY [1..RzextrNmax,1..RzextrNcol] OF real;
   eps,hdid,hnext,htry,x: real;
   i: integer;
   y,dydx,yscal: RealArrayNVAR;

PROCEDURE derivs(x: real; y: RealArrayNVAR; VAR dydx: RealArrayNVAR);
(* Programs using DERIVS must define the type
TYPE
   RealArrayNVAR = 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 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 yscal[i] := 1.0;
   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.