procedure BSSTEP(var Y, DYDX:array of real; NV:integer; var X1:real;
                HTRY, EPS:real;YSCAL:array of real; var HDID, HEND:real);
label 1;
const
    IMAX = 11;  NUSE = 7;  ONE = 1;  SHRINK = 0.95;  GROW = 1.2;
var
    I,J:integer;   H,XSAV,XEST,ERRMAX:real;
    YERR,YSAV,DYSAV,YSEQ:array[0..10] of real;
    NSEQ:array[0..11] of integer;
begin
    NSEQ[1]:=2;   NSEQ[2]:=4;   NSEQ[3]:=6;  NSEQ[4]:=8;  NSEQ[5]:=12;
    NSEQ[6]:=16;  NSEQ[7]:=24;  NSEQ[8]:=32; NSEQ[9]:=48;
    NSEQ[10]:=64; NSEQ[11]:=96;
    H:=HTRY;
    XSAV:=X1;
    For I:=1 To NV do
    begin
        YSAV[I]:=Y[I];
        DYSAV[I]:=DYDX[I];
    end;
1:  For I:=1 To IMAX do
    begin
        MMID(YSAV, DYSAV, NV, XSAV, H, NSEQ[I], YSEQ);
        XEST:=Sqr(H / NSEQ[I]);
        RZEXTR(I, XEST, YSEQ, Y, YERR, NV, NUSE);
        ERRMAX:=0;
        For J:=1 To NV do
            If Abs(YERR[J] / YSCAL[J]) > ERRMAX Then
                ERRMAX:=Abs(YERR[J] / YSCAL[J]);
        ERRMAX:=ERRMAX / EPS;
        If ERRMAX < ONE Then
        begin
            X1:=X1 + H;
            HDID:=H;
            If I = NUSE Then
                HEND:=H * SHRINK
            Else
            begin
                If I = NUSE - 1 Then
                    HEND:=H * GROW
                Else
                    HEND:=(H * NSEQ[NUSE - 1]) / NSEQ[I];
            end;
            Exit;
        end;
    end;
    H:=0.25 * H / Exp(((IMAX - NUSE) div 2)*Ln(2));
    if X1 + H = X1  then  ShowMessage(' Step size underflow');
    goto 1;
end;