procedure ODEINT(var YSTART:array of real; NVAR:integer;X1, X2, EPS,
   H1, HMIN:real;var NOK, NBAD:integer;var XP:array of real;var YP:matrx2);
const
    MAXSTP = 10000;  TWO = 2;  ZERO = 0;  TINY = 1E-30;
var
    YSCAL, Y, DYDX:array[0..10] of real; X,ZZ,H,XSAV,HDID,HEND:real;
    I,NSTP:integer;
begin
    X:=X1;
    If X2-X1 >= 0 THEN
      ZZ:=1
    Else
      ZZ:=-1;
    H:=Abs(H1) * ZZ;
    NOK:=0;
    NBAD:=0;
    KOUNT:=0;
    For I:=1 To NVAR do
        Y[I]:=YSTART[I];
    If KMAX > 0 Then XSAV:=X - DXSAV * TWO;
    For NSTP:=1 To MAXSTP do
    begin
        DERIVS(X, Y, DYDX);
        For I:=1 To NVAR do
            YSCAL[I]:=Abs(Y[I]) + Abs(H * DYDX[I]) + TINY;
        If KMAX > 0 Then
        begin
            If Abs(X - XSAV) > Abs(DXSAV) Then
            begin
                If KOUNT < KMAX - 1 Then
                begin
                    KOUNT:=KOUNT + 1;
                    XP[KOUNT]:=X;
                    For I:=1 To NVAR do
                        YP[I, KOUNT]:=Y[I];
                    XSAV:=X;
                end;
            end;
        end;
        If (X + H - X2) * (X + H - X1) > ZERO Then H:=X2 - X;
        RKQC(Y, DYDX, NVAR, X, H, EPS, YSCAL, HDID, HEND);
        If HDID = H Then 
            NOK:=NOK + 1
        Else
            NBAD:=NBAD + 1;
        If (X - X2) * (X2 - X1) >= ZERO Then
        begin
            For I:=1 To NVAR do
                YSTART[I]:=Y[I];
            If KMAX <> 0 Then
            begin
                KOUNT:=KOUNT + 1;
                XP[KOUNT]:=X;
                For I:=1 To NVAR do
                    YP[I, KOUNT]:=Y[I]
            end;
            Exit;
        end;
        If Abs(HEND) < HMIN Then
        begin
            ShowMessage('Stepsize smaller than minimum.');
            Exit;
        end;
        H:=HEND;
    end; 
    ShowMessage('Too many steps.');
end;