Saturday, November 18, 2017

NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS Y':=F(X,Y) BY ROSENBROCK METHOD.

'****************************************************************
'* NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY *
'* DIFFERENTIAL EQUATIONS Y':=F(X,Y) BY ROSENBROCK METHOD.      *
'* ------------------------------------------------------------ *
'* SAMPLE RUN:                                                  *
'* Example #1:                                                  *
'* (Solve set of differential equations (N=2):                  *
'*  F(1) = Y(1) * Y(2) + COS(X) - HALF * SIN(TWO * X)           *
'*  F(2) = Y(1) * Y(1) + Y(2) * Y(2) - (ONE + SIN(X))           *
'*  Find values of F(1), F(2) at X=1.5).                        *
'*                                                              *
'*  SOLUTION AT X= 1.50000000000000E+0000                       *
'*  Y(1) =   1.236006095804576                                  *
'*  Y(2) =  -.1049268945803322                                  *
'*                                                              *
'*  LAST STEP SIZE = 3.089293500117273D-04                      *
'*  ERROR CODE =  1                                             *
'*                                                              *
'* Example #2:                                                  *
'* (Solve set of differential equations (N=5):                  *
'*  F(1) = Y(2)                                                 *
'*  F(2) = Y(3)                                                 *
'*  F(3) = Y(4)                                                 *
'*  F(4) = Y(5)                                                 *
'*  F(5) = (45.0 * Y(3) * Y(4) * Y(5) -                         *
'*          40.0 * Y(4) * Y(4) * Y(4)) / (9.0 * Y(3) * Y(3))    *
'*  Find values of F(1), F(2), ..., F(5) at X=1.5).             *
'*                                                              *
'*  SOLUTION AT X= 1.50000000000000E+0000                       *
'*  Y(1) =   4.363967162542581                                  *
'*  Y(2) =   4.000019057753676                                  *
'*  Y(3) =   2.82847148934375                                   *
'*  Y(4) =   5.641335228805289D-05                              *
'*  Y(5) =  -3.77130489085171                                   *
'*                                                              *
'*  LAST STEP SIZE =  7.626269049659622D-05                     *
'*  ERROR CODE =  1                                             *
'* ------------------------------------------------------------ *
'* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77   *
'*       [BIBLI 18].                                            *
'*                                                              *
'*                     Basic Release 1.0 By J-P Moreau, Paris   *
'*                               (www.jpmoreau.fr)              *
'****************************************************************
' LIST OF USED SUBROUTINES (HERE INCLUDED):
' ========================================
' 500  FCN      DEFINE SYSTEM OF DIFFERENTIAL EQUATIONS (TWO EXAMPLES)
' 600  IMAX     MAXIMUM OF TWO INTEGERS
' 610  IMIN     MINIMUM OF TWO INTEGERS
' 620  XMAX     MAXIMUM OF TWO REAL NUMBERS
' 630  XMIN     MINIMUM OF TWO REAL NUMBERS
' 650  SIGN     EMULATION OF FUNCTION SIGN OF FORTRAN
' 1000 ROS4     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
'               SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS
'               MY'=F(X,Y)
' 2000 RO4COR   CORE INTEGRATOR FOR ROS4
' 3000 SHAMP    DEFINE CONSTANTS A21... TO D4 (METHOD 1)
' 3001 GRK4A    DEFINE CONSTANTS A21... TO D4 (METHOD 3)
' 3002 GRK4T    DEFINE CONSTANTS A21... TO D4 (METHOD 2 USED HERE)
' 3003 VELDD    DEFINE CONSTANTS A21... TO D4 (METHOD 4)
' 3004 VELDS    DEFINE CONSTANTS A21... TO D4 (METHOD 5)
' 3005 LSTAB    DEFINE CONSTANTS A21... TO D4 (METHOD 6)
' 4000 DECB     MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
'               MATRIX WITH LOWER BANDWIDTH MLDE AND UPPER BANDWIDTH MUE
' 4500 DECA     GENERAL MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION
' 5000 SOLB     SOLUTION OF A BANDED LINEAR SYSTEM, E*X = B
'               (E IS A TRIANGULARIZED MATRIX OBTAINED FROM DECB).
' 5500 SOL      SOLUTION OF A GENERAL LINEAR SYSTEM, E*X = B
'               (E IS A TRIANGULARIZED MATRIX OBTAINED FROM DECA).
' --------------------------------------------------------------------------
'NOTE: THE BANDED MATRIX BRANCH HAS NOT BEEN TESTED HERE, HOWEVER
'      IS FULLY IMPLEMENTED.

DefDbl A-H, O-Z
DefInt I-N

'constants
      NMX = 30  'Maximum size of temporary vectors TMP1 to TMP11 and ITMP

      HALF = 0.5
      ONE = 1#
      TEN = 10#
      TWO = 2#
      XNINE = 9#
      ZERO = 0#

Dim B(NMX), F(NMX), Y(NMX), YY(NMX)
Dim FJAC(5, 5), E(5, 5), FMAS(5, 5)

'variables for statistics (optional use)
'XNFCN,XNSTEP,XNJAC,XNACCPT,NREJCT,XNDEC,XNSOL
'(Long integers are simulated by real numbers).

'begin main program

  'Initialize parameters (see 1000 ROS4)
  N = 2          'DIMENSION OF THE SYSTEM (N=5 for example #2)
  IFCN = 1       'FCN(N,X,Y,F) MAY DEPEND ON X
  X = ZERO       'INITIAL X-VALUE
  XEND = 1.5     'FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
  H = 0.001      'INITIAL STEP SIZE GUESS (0.01 FOR example #2)
  RTOL = 0.00000001 'RELATIVE XERROR TOLERANCE (HERE SCALAR)
  ATOL = 0.0000000001 'ABSOLUTE XERROR TOLERANCE (HERE SCALAR)
  ITOL = 0       'BOTH RTOL AND ATOL ARE SCALARS
  IJAC = 0       'JACOBIAN IS COMPUTED INTERNALLY BY FINITE
                 'DIFFERENCES, Subroutine "JAC" IS NEVER CALLED
  MLJAC = N      'JACOBIAN IS A FULL MATRIX. THE LINEAR ALGEBRA
                 'IS DONE BY FULL-MATRIX GAUSS-ELIMINATION
  IDFX = 0       'DF/DX IS COMPUTED INTERNALLY BY FINITE
                 'DIFFERENCES, Subroutine "DFX" IS NEVER CALLED
  IMAS = 0       'M IS SUPPOSED TO BE THE IDENTITY
                 'MATRIX, Subroutine "MAS" IS NEVER CALLED
  MLMAS = N      'MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ALGEBRA
                 'IS DONE BY FULL-MATRIX GAUSS-ELIMINATION
  IOUT = 0       'Subroutine SOLOUT IS NEVER CALLED
  LE1 = N        'IF MLJAC=N (FULL JACOBIAN)
  LJAC = N       'IF MLJAC=N (FULL JACOBIAN)
  LMAS = 0       'IF IMAS=0

  LIWORK = N + 2                  'DECLARED LENGTH OF ARRAY "IWORK"
  LWORK = N * (LJAC + LMAS + LE1 + 8) + 5 'DECLARED LENGTH OF ARRAY "WORK"

  Dim WORK(LWORK)
  Dim IWORK(LIWORK)

  'Temporary vectors
  Dim TMP1(NMX), TMP2(NMX), TMP3(NMX), TMP4(NMX), TMP5(NMX), TMP6(NMX)
  Dim TMP7(NMX), TMP8(NMX), TMP9(NMX), TMP10(NMX), TMP11(NMX)
  Dim ITMP(NMX)

  For I = 1 To LWORK
    WORK(I) = ZERO 'This triggers default values (see 1000 ROS4)
  Next I

  For I = 1 To LIWORK
    IWORK(I) = 0
  Next I

  Y(1) = HALF    'INITIAL VALUES FOR Y
  Y(2) = HALF    'In example #2, Y(1) = Y(2) = ... = Y(5) = ONE

  Cls
  Print
  Print " Computing..."

  'call Rosenbrock SUBROUTINE with appropriate parameters
  GoSub 1000     'call ROS4(N,IFCN,X,Y,XEND,H,
                 'RTOL,ATOL,ITOL,
                 'IJAC,MLJAC,MUJAC,IDFX,
                 'IMAS,MLMAS,MUMAS,
                 'IOUT,WORK,LWORK,IWORK,LIWORK,IDID)

  'print results
  Cls
  Print
  Print " SOLUTION AT X="; X
  For I = 1 To N
    Print " Y("; I; ") =  "; Y(I)
  Next I
  Print
  Print " LAST STEP SIZE ="; H
  Print " ERROR CODE = "; IDID
  Print
  INPUT "", RR$

End  'of main program


'define example #1
500 'FCN(N,XX,YY,F)
  F(1) = YY(1) * YY(2) + Cos(XX) - HALF * Sin(TWO * XX)
  F(2) = YY(1) * YY(1) + YY(2) * YY(2) - (ONE + Sin(XX))
Return

'define example #2
'500 'FCN(N,XX,YY,F)
'  F(1) = YY(2);
'  F(2) = YY(3);
'  F(3) = YY(4);
'  F(4) = YY(5);
'  F(5) = (45.0 * YY(3) * YY(4) * YY(5) -
'          40.0 * YY(4) * YY(4) * YY(4)) / (NINE * YY(3) * YY(3))
'return


600 'IMAX(ia,ib)
  If ia > ib Then
    IMAX = ia
  Else
    IMAX = ib
  End If
Return

610 'IMIN(ia,ib)
  If ia < ib Then
    IMIN = ia
  Else
    IMIN = ib
  End If
Return

620 'XMAX(xa,xb)
  If xa > xb Then
    XMAX = xa
  Else
    XMAX = xb
  End If
Return

630 'XMIN(xa,xb)
  If xa < xb Then
    XMIN = xa
  Else
    XMIN = xb
  End If
Return

650 'SIGN(xa,xb)
  If xb < 0 Then
    SIGN = -Abs(xa)
  Else
    SIGN = Abs(xa)
  End If
Return



'**********************************************************************
1000 'ROS4 (N,IFCN,X,Y,XEND,H,RTOL,ATOL,ITOL,IJAC,MLJAC,MUJAC,IDFX,
     '          IMAS,MLMAS,MUMAS,IOUT,WORK,LWORK,IWORK,LIWORK,IDID)
' ---------------------------------------------------------------------
'     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
'     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
'     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4
'     (WITH STEP SIZE CONTROL).
'     C.F. SECTION IV.7
'
'     AUTHORS: E. HAIRER AND G. WANNER
'              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
'              CH-1211 GENEVE 24, SWITZERLAND
'              E-MAIL:  HAIRER@CGEUGE51.BITNET,  WANNER@CGEUGE51.BITNET
'
'     THIS CODE IS PART OF THE BOOK:
'         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
'         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
'         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
'         SPRINGER-VERLAG (1990)
'
'     VERSION OF OCTOBER 12, 1990
'
'     INPUT PARAMETERS
'     ----------------
'     N           DIMENSION OF THE SYSTEM
'
'     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
'                 VALUE OF F(X,Y):
'                    Subroutine FCN(N,X,Y,F)
'                    dim Y(N), F(N)
'                    F(1)=...   ETC.
'
'     IFCN        GIVES INFORMATION ON FCN:
'                    IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS)
'                    IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS)
'
'     X           INITIAL X-VALUE
'
'     Y(N)        INITIAL VALUES FOR Y
'
'     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
'
'     H           INITIAL STEP SIZE GUESS;
'                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
'                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
'                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
'                 ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW
'                 STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE.
'                 (IF H=0.D0, THE CODE PUTS H=1.D-6).
'
'     RTOL,ATOL   RELATIVE AND ABSOLUTE XERROR TOLERANCES. THEY
'                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
'
'     ITOL        SWITCH FOR RTOL AND ATOL:
'                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
'                     THE CODE KEEPS, ROUGHLY, THE LOCAL XERROR OF
'                     Y(I) BELOW RTOL*ABS(Y(I))+ATOL
'                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
'                     THE CODE KEEPS THE LOCAL XERROR OF Y(I) BELOW
'                     RTOL(I)*ABS(Y(I))+ATOL(I).
'
'     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
'                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
'                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1).
'                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM:
'                    SUBROUTINE JAC(N,X,Y,DFY,LDFY)
'                    DIM Y(N),DFY(LDFY,N)
'                    DFY(1,1)= ...
'                 LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS
'                 FURNISHED BY THE CALLING PROGRAM.
'                 IF MLJAC = N, THE JACOBIAN IS SUPPOSED TO
'                    BE FULL AND THE PARTIAL DERIVATIVES ARE
'                    STORED IN DFY AS
'                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
'                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
'                    THE PARTIAL DERIVATIVES ARE STORED
'                    DIAGONAL-WISE AS
'                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
'
'     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
'                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
'                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
'                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
'
'     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
'                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
'                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
'                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN
'                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
'                       THE MAIN DIAGONAL).
'
'     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
'                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
'                 NEED NOT BE DEFINED IF MLJAC=N.
'
'     DFX         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
'                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X
'                 (THIS ROUTINE IS ONLY CALLED IF IDFX=1 AND IFCN=1).
'                 OTHERWISE, THIS SUBROUTINE MUST HAVE THE FORM
'                    SUBROUTINE DFX(N,X,Y,FX)
'                    DIM Y(N),FX(N)
'                    FX(1)= ...
'
'     IDFX        SWITCH FOR THE COMPUTATION OF THE DF/DX:
'                    IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE
'                       DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED.
'                    IDFX=1: DF/DX IS SUPPLIED BY SUBROUTINE DFX.
'
'     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
'     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
'
'     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
'                 MATRIX M.
'                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
'                 MATRIX AND NEEDS NOT TO BE DEFINED;
'                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
'                    SUBROUTINE MAS(N,AM,LMAS)
'                    DIM AM(LMAS,N)
'                    AM(1,1)= ....
'                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
'                    AS FULL MATRIX LIKE
'                         AM(I,J) = M(I,J)
'                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
'                    DIAGONAL-WISE AS
'                         AM(I-J+MUMAS+1,J) = M(I,J).
'
'     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
'                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
'                       MATRIX, MAS IS NEVER CALLED.
'                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
'
'     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
'                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
'                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
'                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
'                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
'                       THE MAIN DIAGONAL).
'                 MLMAS IS SUPPOSED TO BE <= MLJAC.
'
'     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
'                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
'                 NEED NOT BE DEFINED IF MLMAS=N.
'                 MUMAS IS SUPPOSED TO BE <= MUJAC.
'
'     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
'                 NUMERICAL SOLUTION DURING INTEGRATION.
'                 IF IOUT=0, NO INTERPOLATION SUBROUTINE IS NECESSARY.
'                 IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
'                 IT MUST HAVE THE FORM
'                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN)
'                    DIM Y(N)
'                    ....
'                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
'                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
'                    THE FIRST GRID-POINT).
'                 "IRTRN" SERVES TO INTXERRUPT THE INTEGRATION. IF IRTRN
'                    IS SET <0, ROS4 RETURNS TO THE CALLING PROGRAM.
'
'     IOUT        GIVES INFORMATION ON THE SUBROUTINE SOLOUT:
'                    IOUT=0: SUBROUTINE IS NEVER CALLED
'                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
'
'     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
'                 SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES.
'                 "LWORK" MUST BE AT LEAST
'                             N*(LJAC+LMAS+LE1+8)+5
'                 WHERE
'                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
'                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
'                 AND
'                    LMAS=0              IF IMAS=0
'                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
'                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
'                 AND
'                    LE1=N               IF MLJAC=N (FULL JACOBIAN)
'                    LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.).
'
'                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
'                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
'                 STORAGE REQUIREMENT IS
'                             LWORK = 2*N*N+8*N+5.
'
'     LWORK       DECLARED LENGHT OF ARRAY "WORK".
'
'     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
'                 "LIWORK" MUST BE AT LEAST N+2.
'
'     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
'
' ----------------------------------------------------------------------
'
'     SOPHISTICATED SETTING OF PARAMETERS
'     -----------------------------------
'              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
'              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(5)
'              AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO.
'              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
'
'    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
'              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
'
'    IWORK(2)  SWITCH FOR THE CHOICE OF THE COEFFICIENTS
'              IF IWORK(2).EQ.1  METHOD OF SHAMPINE
'              IF IWORK(2).EQ.2  METHOD GRK4T OF KAPS-RENTROP
'              IF IWORK(2).EQ.3  METHOD GRK4A OF KAPS-RENTROP
'              IF IWORK(2).EQ.4  METHOD OF VAN VELDHUIZEN (GAMMA=1/2)
'              IF IWORK(2).EQ.5  METHOD OF VAN VELDHUIZEN ("D-STABLE")
'              IF IWORK(2).EQ.6  AN L-STABLE METHOD
'              THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=2.
'
'    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1D-16.
'
'    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
'
'    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
'              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
'                 WORK(3) <= HNEW/HOLD <= WORK(4)
'              DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0
'
'    WORK(5)   AVOID THE HUMP: AFTER TWO CONSECUTIVE STEP REJECTIONS
'              THE STEP SIZE IS MULTIPLIED BY WORK(5)
'              DEFAULT VALUES: WORK(5)=0.1#
'
'-----------------------------------------------------------------------
'
'     OUTPUT PARAMETERS
'     -----------------
'     X           X-VALUE WHERE THE SOLUTION IS COMPUTED
'                 (AFTER SUCCESSFUL RETURN X=XEND)
'
'     Y(N)       SOLUTION AT X
'
'     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
'
'     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
'                   IDID=1  COMPUTATION SUCCESSFUL,
'                   IDID=-1 COMPUTATION UNSUCCESSFUL.
'
' ---------------------------------------------------------
' *** *** *** *** *** *** *** *** *** *** *** *** ***
'          DECLARATIONS
' *** *** *** *** *** *** *** *** *** *** *** *** ***

' IAUTNMS,IMPLCT,JBAND,IARRET: Boolean (here integers).
' --------------------------------------------------------------------
' --- THESE COMMON VARIABLES CAN BE USED FOR STATISTICS:
' ---    XNFCN     NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
'                  EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
' ---    XNJAC     NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
'                  OR NUMERICALLY)
' ---    XNSTEP    NUMBER OF COMPUTED STEPS
' ---    XNACCPT   NUMBER OF ACCEPTED STEPS
' ---    NREJCT    NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
'                  HAS BEEN ACCEPTED)
' ---    XNDEC     NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
' ---    XNSOL     NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
' --------------------------------------------------------------------
'     LDE,LDJAC,LDMAS,LDMAS2,METH: Integer
'     XNMAX: LongInt (here double)
'     FAC1,FAC2,FACREJ,HMAX,UROUND: Double
'     IEYNEW,IEDY1,IEDY,IEAK1,IEAK2,IEAK3,IEAK4,
'     IEFX,IEJAC,IEMAS, IEE, ISTORE: Integer;
'     IEIP: Integer

' *** *** *** *** *** *** ***
'    SETTING THE PARAMETERS
' *** *** *** *** *** *** ***
      XNFCN = ZERO
      XNJAC = ZERO
      XNSTEP = ZERO
      XNACCPT = ZERO
      NREJCT = 0
      XNDEC = ZERO
      XNSOL = ZERO
      IARRET = 0
' ------- XNMAX , THE MAXIMAL NUMBER OF STEPS -----
      If IWORK(1) = 0 Then
         XNMAX = 100000#
      Else
         XNMAX = 1# * IWORK(1)
         If XNMAX <= ZERO Then
           Print " WRONG INPUT, IWORK(1)= "; IWORK(1)
           IARRET = 1
         End If
      End If
' -------- METH   COEFFICIENTS OF THE METHOD ------
      If IWORK(2) = 0 Then
         METH = 2
      Else
         METH = IWORK(2)
         If METH <= 0 Or METH >= 7 Then
           Print " CURIOUS INPUT, IWORK(2)="; IWORK(2)
           IARRET = 1
         End If
      End If
' -------- UROUND, SMALLEST NUMBER SATISFYING ONE + UROUND > ONE ---
      If WORK(1) = ZERO Then
         UROUND = 1E-16
      Else
         UROUND = WORK(1)
         If UROUND <= 0.00000000000001 Or UROUND >= ONE Then
           Print " COEFFICIENTS HAVE 16 DIGITS, UROUND="; WORK(1)
           IARRET = 1
         End If
      End If
' -------- MAXIMAL STEP SIZE -----------
      If WORK(2) = ZERO Then
         HMAX = XEND - X
      Else
         HMAX = WORK(2)
      End If
' -------  FAC1,FAC2   PARAMETERS FOR STEP SIZE SELECTION ------
      If WORK(3) = ZERO Then
         FAC1 = 5#
      Else
         FAC1 = ONE / WORK(3)
      End If
      If WORK(4) = ZERO Then
         FAC2 = ONE / 6#
      Else
         FAC2 = ONE / WORK(4)
      End If
' -------  FACREJ    FOR THE HUMP -------
      If WORK(5) = ZERO Then
         FACREJ = 0.1
      Else
         FACREJ = WORK(5)
      End If
' --------- CHECK IF TOLERANCES ARE O.K. ---------
      If ITOL = 0 Then
          If ATOL <= ZERO Or RTOL <= TEN * UROUND Then
            Print " TOLERANCES ARE TOO SMALL."
            IARRET = 1
          End If
      Else
         'Multiple tolerances not implemented here
         'For I=1 to N
         'IF ATOL(I) <= ZERO OR RTOL(I) <= TEN*UROUND THEN
         '  print " TOLERANCE(";I;") IS TOO SMALL."
         '  IARRET=1
         'end if
      End If
' *** *** *** *** *** *** *** *** *** *** *** *** ***
'         COMPUTATION OF ARRAY ENTRIES
' *** *** *** *** *** *** *** *** *** *** *** *** ***
' ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
      If IFCN = 0 Then
        IAUTNMS = 1
      Else
        IAUTNMS = 0
      End If
      If IMAS <> 0 Then
        IMPLCT = 1
      Else
        IMPLCT = 0
      End If
      If MLJAC <> N Then
        JBAND = 1
      Else
        JBAND = 0
      End If

      IARRET = 0

' -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ------
' -- JACOBIAN --
      If JBAND <> 0 Then
         LDJAC = MLJAC + MUJAC + 1
      Else
         LDJAC = N
      End If
' -- MATRIX E FOR LINEAR ALGEBRA --
      If JBAND <> 0 Then
         LDE = 2 * MLJAC + MUJAC + 1
      Else
         LDE = N
      End If
' -- MASS MATRIX --
      If IMPLCT <> 0 Then
        If MLMAS <> N Then
          LDMAS = MLMAS + MUMAS + 1
        Else
          LDMAS = N
        End If
' ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC" -------
        If MLMAS > MLJAC Or MUMAS > MUJAC Then
          Print " BANDWITH OF MAS MUST NOT BE LARGER THAN BANDWITH OF JAC."
          IARRET = 1
        End If
      Else
        LDMAS = 0
      End If

      ia = 1: ib = LDMAS: GoSub 600
      LDMAS2 = IMAX

' ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ------
      IEYNEW = 6
      IEDY1 = IEYNEW + N
      IEDY = IEDY1 + N
      IEAK1 = IEDY + N
      IEAK2 = IEAK1 + N
      IEAK3 = IEAK2 + N
      IEAK4 = IEAK3 + N
      IEFX = IEAK4 + N
      IEJAC = IEFX + N
      IEMAS = IEJAC + N * LDJAC
      IEE = IEMAS + N * LDMAS
' ------ TOTAL STORAGE REQUIREMENT -----------
      ISTORE = IEE + N * LDE - 1
      If ISTORE > LWORK Then
        Print " INSUFFICIENT STORAGE FOR WORK, MINIMUM LWORK="; ISTORE
        IARRET = 1
      End If
' ------- ENTRY POINTS FOR INTEGER WORKSPACE ------
      IEIP = 3
' ---------  TOTAL REQUIREMENT  ---------------
      ISTORE = IEIP + N - 1
      If ISTORE > LIWORK Then
        Print " INSUFF. STORAGE FOR IWORK, MINIMUM LIWORK="; ISTORE
        IARRET = 1
      End If
' ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 ------
      If IARRET <> 0 Then
         IDID = -1
         Return
      End If
' ------------- prepare arguments of RO4COR --------------------------
' Here, appropriate parts of WORK (IWORK) are put in temporary vectors
' TMP1 to TMP11 (ITMP) to simulate Fortran arguments of RO4COR.
      For I = 1 To N
        If I + IEYNEW - 1 <= N Then
          TMP1(I) = WORK(I + IEYNEW - 1)
        Else
          TMP1(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEDY1 - 1 <= N Then
          TMP2(I) = WORK(I + IEDY1 - 1)
        Else
          TMP2(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEDY - 1 <= N Then
          TMP3(I) = WORK(I + IEDY - 1)
        Else
          TMP3(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEAK1 - 1 <= N Then
          TMP4(I) = WORK(I + IEAK1 - 1)
        Else
          TMP4(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEAK2 - 1 <= N Then
          TMP5(I) = WORK(I + IEAK2 - 1)
        Else
          TMP5(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEAK3 - 1 <= N Then
          TMP6(I) = WORK(I + IEAK3 - 1)
        Else
          TMP6(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEAK4 - 1 <= N Then
          TMP7(I) = WORK(I + IEAK4 - 1)
        Else
          TMP7(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEFX - 1 <= N Then
          TMP8(I) = WORK(I + IEFX - 1)
        Else
          TMP8(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEJAC - 1 <= N Then
          TMP9(I) = WORK(I + IEJAC - 1)
        Else
          TMP9(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEE - 1 <= N Then
          TMP10(I) = WORK(I + IEE - 1)
        Else
          TMP10(I) = ZERO
        End If
      Next I
      For I = 1 To No
        If I + IEMAS - 1 <= N Then
          TMP11(I) = WORK(I + IEMAS - 1)
        Else
          TMP11(I) = ZERO
        End If
      Next I
      For I = 1 To N
        If I + IEIP - 1 <= N Then
          ITMP(I) = IWORK(I + IEIP - 1)
        Else
          ITMP(I) = 0
        End If
      Next I

' --------  CALL TO CORE INTEGRATOR ------------
      GoSub 2000 'call RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IJAC,
                 'MLJAC,MUJAC,IDFX,MLMAS,MUMAS,IOUT,IDID,
                 'XNMAX,UROUND,METH,FAC1,FAC2,FACREJ,IAUTNMS,
                 'IMPLCT,JBAND,LDJAC,LDE,LDMAS2,TMP1,TMP2,
                 'TMP3,TMP4,TMP5,TMP6,TMP7,TMP8,TMP9,TMP10,
                 'TMP11, ITMP)

Return 'from 1000 ROS4


'  --------- ...  AND HERE IS THE CORE INTEGRATOR  ----------

2000 'Subroutine RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IJAC,
                'MLJAC,MUJAC,IDFX,MLMAS,MUMAS,IOUT,IDID,
                'XNMAX,UROUND,METH,FAC1,FAC2,FACREJ,IAUTNMS,
                'IMPLCT,JBAND,LDJAC,LDE,LDMAS2,TMP1,TMP2,
                'TMP3,TMP4,TMP5,TMP6,TMP7,TMP8, TMP9,TMP10,
                'TMP11,ITMP)
' ----------------------------------------------------------
'     CORE INTEGRATOR FOR ROS4
'     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED
' ----------------------------------------------------------
'         DECLARATIONS
' ----------------------------------------------------------
'Labels 2001,2002, 2012, 2014, 2079, 2080

'     I,J,K,L,MBB,MBDIAG,MBJAC,MDIAG,MDIFF,MLDE,MUE: Integer
'     IREJECT,IRJECT2: Boolean (here integers)
'     A21,A31,A32,C21,C31,C32,C41,C42,C43: Double
'     B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4: Double
'     DELT,HMAXN,HOPT,POSNEG,XDELT,XXOLD,YSAFE: Double
'     IRTRN,LBEG,LDEND,MD,MUJACJ,MUJACP,NSING:Integer;
'     FAC,HC21,HC31,HC32,HC41,HC42,HC43,SUM: Double
'     I1,I2,IB,INFO,J1,J2,MADD: Integer;
'     XERR, HD1,HD2,HD3,HD4, HNEW,S,SK: Double

' ------- restore Fortran parameters FJAC, E, FMAS ----------
      For J = 1 To N
        For I = 1 To N
          FJAC(I, J) = TMP9(I + J - 1)
          E(I, J) = TMP10(I + J - 1)
          FMAS(I, J) = TMP11(I + J - 1)
        Next I
      Next J
' ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
'     IF IMPLCT <> 0 CALL MAS(N,FMAS,LDMAS2)
'     (Not provided here).
' ---- PREPARE BANDWIDTHS -----
      If JBAND <> 0 Then
        MLDE = MLJAC
        MUE = MUJAC
        MBJAC = MLJAC + MUJAC + 1
        MBB = MLMAS + MUMAS + 1
        MDIAG = MLDE + MUE + 1
        MBDIAG = MUMAS + 1
        MDIFF = MLDE + MUE - MUMAS
      End If
' *** *** *** *** *** *** ***
'  INITIALISATIONS
' *** *** *** *** *** *** *** }
      'POSNEG=SIGN(ONE,XEND-X)
      xa = ONE: xb = XEND - X: GoSub 650: POSNEG = SIGN

      If METH = 1 Then GoSub 3000 'call SHAMP(A21,A31,A32,C21,C31,C32,C41,C42,C43,
                                  'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      If METH = 2 Then GoSub 3002 'call GRK4T(A21,A31,A32,C21,C31,C32,C41,C42,C43,
                                  'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      If METH = 3 Then GoSub 3001 'call GRK4A(A21,A31,A32,C21,C31,C32,C41,C42,C43,
                                  'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      If METH = 4 Then GoSub 3003 'call VELDS(A21,A31,A32,C21,C31,C32,C41,C42,C43,
                                  'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      If METH = 5 Then GoSub 3004 'call VELDD(A21,A31,A32,C21,C31,C32,C41,C42,C43,
                                  'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      If METH = 6 Then GoSub 3005 'call LSTAB(A21,A31,A32,C21,C31,C32,C41,C42,C43,
                                  'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)

' --- INITIAL PREPARATIONS ---

      'HMAXN=XMIN(ABS(HMAX),ABS(XEND-X))
      xa = Abs(HMAX): xb = Abs(XEND - X): GoSub 630
      HMAXN = XMIN
      'H=XMIN(XMAX(1D-10,ABS(H)),HMAXN)
      xa = 0.0000000001: xb = Abs(H): GoSub 620
      xa = XMAX: xb = HMAXN: GoSub 630: H = XMIN
      'H=SIGN(H,POSNEG)
      xa = H: xb = POSNEG: GoSub 650: H = SIGN
      IREJECT = 0
      NSING = 0
      IRTRN = 1
      XXOLD = X

'     IF IOUT <> 0 THEN
'       call SOLOUT(XNACCPT+1,XXOLD,X,Y,N,IRTRN)
'       (not provided here).
'     END IF
     
      If IRTRN < 0 Then GoTo 2079

' --- BASIC INTEGRATION STEP ---

2001  If XNSTEP > XNMAX Or X + 0.1 * H = X Or Abs(H) <= UROUND Then GoTo 2079
      If (X - XEND) * POSNEG + UROUND > ZERO Then
        H = HOPT
        IDID = 1
        Return    'Normal return
      End If
      HOPT = H
      If (X + H - XEND) * POSNEG > ZERO Then H = XEND - X

      'call FCN(N,X,Y,TMP2)
      XX = X
      For I = 1 To N
        YY(I) = Y(I)
      Next I
      GoSub 500
      For I = 1 To N
        TMP2(I) = F(I)
      Next I

      XNFCN = XNFCN + 1#
'  *** *** *** *** *** *** ***
'  COMPUTATION OF THE JACOBIAN
'  *** *** *** *** *** *** ***
      XNJAC = XNJAC + ONE
      If IJAC = 0 Then
' --- COMPUTE JACOBIAN MATRIX NUMERICALLY ---
          If JBAND <> 0 Then
' --- JACOBIAN IS BANDED ---
            MUJACP = MUJAC + 1
            'MD=IMIN(MBJAC,N)
            ia = MBJAC: ib = N: GoSub 610: MD = IMIN
            For K = 1 To MD
              J = K
2012          TMP5(J) = Y(J)
              'TMP6(J)=SQR(UROUND*XMAX(1D-5,ABS(Y(J))))
              xa = 0.00001: xb = Y(J): GoSub 620
              TMP6(J) = Sqr(UROUND * XMAX)
              Y(J) = Y(J) + TMP6(J)
              J = J + MD
              If J <= N Then GoTo 2012
              'call FCN(N,X,Y,TMP4)
              GoSub 500
              For I = 1 To N
                TMP4(I) = F(I)
              Next I
              J = K
              'LBEG=IMAX(1,J-MUJAC)
              ia = 1: ib = MUJAC: GoSub 600: LBEG = IMAX
2014          'LDEND=IMIN(N,J+MLJAC)
              ia = N: ib = J + MLJAC: GoSub 610: LDEND = IMIN
              Y(J) = TMP5(J)
              MUJACJ = MUJACP - J
              For L = LBEG To LDEND
                FJAC(L + MUJACJ, J) = (TMP4(L) - TMP2(L)) / TMP6(J)
              Next L
              J = J + MD
              LBEG = LDEND + 1
              If J <= N Then GoTo 2014
            Next K
          Else
' --- JACOBIAN IS FULL ---
            For I = 1 To N
              YSAFE = Y(I)
              'DELT=SQR(UROUND*XMAX(1D-5,ABS(YSAFE)))
              xa = 0.00001: xb = Abs(YSAFE): GoSub 620
              DELT = Sqr(UROUND * XMAX)
              Y(I) = YSAFE + DELT
              'call FCN(N,X,Y,TMP4)
              XX = X
              For II = 1 To N
                YY(II) = Y(II)
              Next II
              GoSub 500
              For II = 1 To N
                TMP4(II) = F(II)
              Next II
              For J = 1 To N
                FJAC(J, I) = (TMP4(J) - TMP2(J)) / DELT
              Next J
              Y(I) = YSAFE
            Next I
            MLJAC = N
          End If
      Else
' --- COMPUTE JACOBIAN MATRIX ANALYTICALLY ---
'         JAC(N,X,Y,FJAC,LDJAC)
'         (Not provided here).
      End If
      If IAUTNMS = 0 Then
          If IDTMP8 = 0 Then
' --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X ---
            'DELT=SQR(UROUND*XMAX(1D-5,ABS(X)))
            xa = 0.00001: xb = Abs(X): GoSub 620
            DELT = Sqr(UROUND * XMAX)
            XDELT = X + DELT
            'call FCN(N,XDELT,Y,TMP4)
            XX = XDELT: GoSub 500
            For I = 1 To N
              TMP4(I) = F(I)
            Next I
            For J = 1 To N
              TMP8(J) = (TMP4(J) - TMP2(J)) / DELT
            Next J
'         ELSE
' --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X ---
'           CALL DTMP8(N,X,Y,TMP8)
'           (Not provided here).
          End If
      End If
' *** *** *** *** *** *** ***
'      COMPUTE THE STAGES
' *** *** *** *** *** *** ***
2002  XNDEC = XNDEC + ONE
      HC21 = C21 / H
      HC31 = C31 / H
      HC32 = C32 / H
      HC41 = C41 / H
      HC42 = C42 / H
      HC43 = C43 / H
      FAC = ONE / (H * GAMMA)
      If IMPLCT <> 0 Then
          If JBAND <> 0 Then
' --- THE MATRIX E (B IS A JBAND MATRIX, JACOBIAN A JBAND MATRIX) ---
            For J = 1 To N
              'I1=IMAX(1,MUJAC+2-J)
              ia = 1: ib = MUJAC + 2 - J: GoSub 600: I1 = IMAX
              'I2=IMIN(MBJAC,N+MUJAC+1-J)
              ia = MBJAC: ib = N + MUJAC + 1 - J: GoSub 610: I2 = IMIN
              For I = I1 To I2
                E(I + MLDE, J) = -FJAC(I, J)
              Next I
            Next J
            For J = 1 To N
              'I1=IMAX(1,MUMAS+2-J)
              ia = 1: ib = MUMAS + 2 - J: GoSub 600: I1 = IMAX
              'I2=IMIN(MBB,N+MUMAS+1-J)
              ia = MBB: ib = N + MUMAS + 1 - J: GoSub 610: I2 = IMIN
              For I = I1 To I2
                ib = I + MDIFF
                E(ib, J) = E(ib, J) + FAC * FMAS(I, J)
              Next I
            Next J
            GoSub 4000  'call DECB(N,E,MLDE,MUE,ITMP,INFO)
            If INFO <> 0 Then GoTo 2080
            If IAUTNMS <> 0 Then
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
' ---   2) THE MATRIX B AND THE JACOBIAN OF F ARE BANDED
' ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS.
                For I = 1 To N
                  TMP4(I) = TMP2(I)
                Next I
                For I = 1 To N
                  B(I) = TMP4(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP)
                For I = 1 To N
                  TMP4(I) = B(I)
                Next I
                For I = 1 To N
                  TMP1(I) = Y(I) + A21 * TMP4(I)
                Next I
                'call FCN(N,X,TMP1,TMP3)
                XX = X
                For I = 1 To N
                  YY(I) = TMP1(I)
                Next I
                GoSub 500
                For I = 1 To N
                  TMP3(I) = F(I)
                Next I
                For I = 1 To N
                  TMP1(I) = HC21 * TMP4(I)
                Next I
                For I = 1 To N
                  Sum = ZERO
                  'J1=IMAX(1,I-MLMAS)
                  ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                  'J2=IMIN(N,I+MUMAS)
                  ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                  For J = J1 To J2
                    Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J)
                  Next J
                  TMP5(I) = Sum + TMP3(I)
                Next I
                For I = 1 To N
                  B(I) = TMP5(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP)
                For I = 1 To N
                  TMP5(I) = B(I)
                Next I
                For I = 1 To N
                  TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
                Next I
                'call FCN(N,X,TMP1,TMP3)
                For I = 1 To N
                  YY(I) = TMP1(I)
                Next I
                GoSub 500
                For I = 1 To N
                  TMP3(I) = F(I)
                Next I
                For I = 1 To N
                  TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I)
                Next I
                For I = 1 To N
                  Sum = ZERO
                  'J1=IMAX(1,I-MLMAS)
                  ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                  'J2=IMIN(N,I+MUMAS)
                  ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                  For J = J1 To J2
                    Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J)
                  Next J
                  TMP6(I) = Sum + TMP3(I)
                Next I
                For I = 1 To N
                  B(I) = TMP6(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP)
                For I = 1 To N
                  TMP6(I) = B(I)
                Next I
                For I = 1 To N
                  TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
                Next I
                For I = 1 To N
                  Sum = ZERO
                  'J1=IMAX(1,I-MLMAS)
                  ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                  'J2=IMIN(N,I+MUMAS)
                  ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                  For J = J1 To J2
                    Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J)
                  Next J
                  TMP7(I) = Sum + TMP3(I)
                Next I
                For I = 1 To N
                  B(I) = TMP7(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP)
                For I = 1 To N
                  TMP7(I) = B(I)
                Next I
            Else
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
' ---   2) THE MATRIX B AND THE JACOBIAN OF F ARE BANDED
' ---   3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS.
                HD1 = H * D1
                HD2 = H * D2
                HD3 = H * D3
                HD4 = H * D4
                For I = 1 To N
                  TMP4(I) = TMP2(I) + HD1 * TMP8(I)
                Next I
                For I = 1 To N
                  B(I) = TMP4(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP)
                For I = 1 To N
                  TMP4(I) = B(I)
                Next I
                For I = 1 To N
                  TMP1(I) = Y(I) + A21 * TMP4(I)
                Next I
                'call FCN(N,X+C2*H,TMP1,TMP3)
                XX = X + C2 * H
                For I = 1 To N
                  YY(I) = TMP1(I)
                Next I
                GoSub 500
                For I = 1 To N
                  TMP3(I) = F(I)
                Next I
                For I = 1 To N
                  TMP1(I) = HC21 * TMP4(I)
                Next I
                For I = 1 To N
                  Sum = ZERO
                  'J1=IMAX(1,I-MLMAS)
                  ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                  'J2=IMIN(N,I+MUMAS)
                  ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                  For J = J1 To J2
                    Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J)
                  Next J
                  TMP5(I) = Sum + TMP3(I) + HD2 * TMP8(I)
                Next I
                For I = 1 To N
                  B(I) = TMP5(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP)
                For I = 1 To N
                  TMP5(I) = B(I)
                Next I
                For I = 1 To N
                  TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
                Next I
                'call FCN(N,X+C3*H,TMP1,TMP3)
                XX = X + C3 * H
                For I = 1 To N
                  YY(I) = TMP1(I)
                Next I
                GoSub 500
                For I = 1 To N
                  TMP3(I) = F(I)
                Next I
                For I = 1 To N
                  TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I)
                Next I
                For I = 1 To N
                  Sum = ZERO
                  'J1=IMAX(1,I-MLMAS)
                  ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                  'J2=IMIN(N,I+MUMAS)
                  ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                  For J = J1 To J2
                    Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J)
                  Next J
                  TMP6(I) = Sum + TMP3(I) + HD3 * TMP8(I)
                Next I
                For I = 1 To N
                  B(I) = TMP6(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP)
                For I = 1 To N
                  TMP6(I) = B(I)
                Next I
                For I = 1 To N
                  TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
                Next I
                For I = 1 To N
                  Sum = ZERO
                  'J1=IMAX(1,I-MLMAS)
                  ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                  'J2=IMIN(N,I+MUMAS)
                  ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                  For J = J1 To J2
                    Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J)
                  Next J
                  TMP7(I) = Sum + TMP3(I) + HD4 * TMP8(I)
                Next I
                For I = 1 To N
                  B(I) = TMP7(I)
                Next I
                GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP)
                For I = 1 To N
                  TMP7(I) = B(I)
                Next I
            End If
          Else
              If MLMAS <> N Then
' --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX) ---
                MADD = MUMAS + 1
                For J = 1 To N
                  For I = 1 To N
                    E(I, J) = -FJAC(I, J)
                  Next I
                Next J
                For J = 1 To N
                  'I1=IMAX(1,J-MUMAS)
                  ia = 1: ib = J - MUMAS: GoSub 600: I1 = IMAX
                  'I2=IMIN(N,J+MLMAS)
                  ia = N: ib = J + MLMAS: GoSub 610: I2 = IMIN
                  For I = I1 To I2
                    E(I, J) = E(I, J) + FAC * FMAS(I - J + MADD, J)
                  Next I
                Next J
                GoSub 4500 'call DECA(N,LDE,E,ITMP,INFO)
                If INFO <> 0 Then GoTo 2080
                If IAUTNMS <> 0 Then
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
' ---   2) THE MATRIX B IS BANDED BUT THE JACOBIAN OF F IS NOT
' ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS.
                    For I = 1 To N
                      TMP4(I) = TMP2(I)
                    Next I
                    For I = 1 To N
                      B(I) = TMP4(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP4,ITMP)
                    For I = 1 To N
                      TMP4(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A21 * TMP4(I)
                    Next I
                    'call FCN(N,X,TMP1,TMP3)
                    XX = X
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC21 * TMP4(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I)
                      'J1=IMAX(1,I-MLMAS)
                      ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                      'J2=IMIN(N,I+MUMAS)
                      ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                      For J = J1 To J2
                        Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J)
                      Next J
                      TMP5(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP5(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP5,ITMP)
                    For I = 1 To N
                      TMP5(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
                    Next I
                    'call FCN(N,X,TMP1,TMP3)
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I)
                      'J1=IMAX(1,I-MLMAS)
                      ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                      'J2=IMIN(N,I+MUMAS)
                      ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                      For J = J1 To J2
                        Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J)
                      Next J
                      TMP6(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP6(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP6,ITMP)
                    For I = 1 To N
                      TMP6(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I)
                      'J1=IMAX(1,I-MLMAS)
                      ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                      'J2=IMIN(N,I+MUMAS)
                      ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                      For J = J1 To J2
                        Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J)
                      Next J
                      TMP7(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP7(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP7,ITMP)
                    For I = 1 To N
                      TMP7(I) = B(I)
                    Next I
                Else
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
' ---   2) THE MATRIX B IS BANDED BUT THE JACOBIAN OF F IS NOT
' ---   3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS.
                    HD1 = H * D1
                    HD2 = H * D2
                    HD3 = H * D3
                    HD4 = H * D4
                    For I = 1 To N
                      TMP4(I) = TMP2(I) + HD1 * TMP8(I)
                    Next I
                    For I = 1 To N
                      B(I) = TMP4(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP4,ITMP)
                    For I = 1 To N
                      TMP4(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A21 * TMP4(I)
                    Next I
                    'call FCN(N,X+C2*H,TMP1,TMP3)
                    XX = X + C2 * H
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC21 * TMP4(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I) + HD2 * TMP8(I)
                      'J1=IMAX(1,I-MLMAS)
                      ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                      'J2=IMIN(N,I+MUMAS)
                      ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                      For J = J1 To J2
                        Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J)
                      Next J
                      TMP5(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP5(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP5,ITMP)
                    For I = 1 To N
                      TMP5(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
                    Next I
                    'call FCN(N,X+C3*H,TMP1,TMP3)
                    XX = X + C3 * H
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I) + HD3 * TMP8(I)
                      'J1=IMAX(1,I-MLMAS)
                      ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                      'J2=IMIN(N,I+MUMAS)
                      ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                      For J = J1 To J2
                        Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J)
                      Next J
                      TMP6(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP5(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP5,ITMP)
                    For I = 1 To N
                      TMP5(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I) + HD4 * TMP8(I)
                      'J1=IMAX(1,I-MLMAS)
                      ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX
                      'J2=IMIN(N,I+MUMAS)
                      ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN
                      For J = J1 To J2
                        Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J)
                      Next J
                      TMP7(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP7(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP7,ITMP)
                    For I = 1 To N
                      TMP7(I) = B(I)
                    Next I
                End If
              Else
' --- THE MATRIX E (B IS A FULL MATRIX, JACOBIAN A FULL OR BANDED MATRIX) ---
                  If MLJAC = N Then
                    For J = 1 To N
                      For I = 1 To N
                        E(I, J) = FMAS(I, J) * FAC - FJAC(I, J)
                      Next I
                    Next J
                  Else
                    MADD = MUJAC + 1
                    For J = 1 To N
                      For I = 1 To N
                        E(I, J) = FMAS(I, J) * FAC
                      Next I
                    Next J
                    For J = 1 To N
                      'I1=IMAX(1,J-MUJAC)
                      ia = 1: ib = J - MUJAC: GoSub 600: I1 = IMAX
                      'I2=IMIN(N,J+MLJAC)
                      ia = N: ib = J + MLJAC: GoSub 610: I2 = IMIN
                      For I = I1 To I2
                        E(I, J) = E(I, J) - FJAC(I - J + MADD, J)
                      Next I
                    Next J
                  End If
                  GoSub 4500 'call DECA(N,E,ITMP,INFO)
                  If INFO <> 0 Then GoTo 2080
                  If IAUTNMS <> 0 Then
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
' ---   2) THE MATRIX B IS NOT BANDED
' ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS.
                    For I = 1 To N
                      TMP4(I) = TMP2(I)
                    Next I
                    For I = 1 To N
                      B(I) = TMP4(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP4,ITMP)
                    For I = 1 To N
                      TMP4(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A21 * TMP4(I)
                    Next I
                    'call FCN(N,X,TMP1,TMP3)
                    XX = X
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC21 * TMP4(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I)
                      For J = 1 To N
                        Sum = Sum + FMAS(I, J) * TMP1(J)
                      Next J
                      TMP5(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP5(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP5,ITMP)
                    For I = 1 To N
                      TMP5(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
                    Next I
                    'call FCN(N,X,TMP1,TMP3)
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I)
                      For J = 1 To N
                        Sum = Sum + FMAS(I, J) * TMP1(J)
                      Next J
                      TMP6(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP6(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP6,ITMP)
                    For I = 1 To N
                      TMP6(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I)
                      For J = 1 To N
                        Sum = Sum + FMAS(I, J) * TMP1(J)
                      Next J
                      TMP7(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP7(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP7,ITMP)
                    For I = 1 To N
                      TMP7(I) = B(I)
                    Next I
                  Else
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
' ---   2) THE MATRIX B IS NOT BANDED
' ---   3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS.
                    HD1 = H * D1
                    HD2 = H * D2
                    HD3 = H * D3
                    HD4 = H * D4
                    For I = 1 To N
                      TMP4(I) = TMP2(I) + HD1 * TMP8(I)
                    Next I
                    For I = 1 To N
                      B(I) = TMP4(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP4,ITMP)
                    For I = 1 To N
                      TMP4(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A21 * TMP4(I)
                    Next I
                    'call FCN(N,X+C2*H,TMP1,TMP3)
                    XX = X + C2 * H
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC21 * TMP4(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I) + HD2 * TMP8(I)
                      For J = 1 To N
                        Sum = Sum + FMAS(I, J) * TMP1(J)
                      Next J
                      TMP5(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP5(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP5,ITMP)
                    For I = 1 To N
                      TMP5(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
                    Next I
                    'call FCN(N,X+C3*H,TMP1,TMP3)
                    XX = X + C3 * H
                    For I = 1 To N
                      YY(I) = TMP1(I)
                    Next I
                    GoSub 500
                    For I = 1 To N
                      TMP3(I) = F(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I) + HD3 * TMP8(I)
                      For J = 1 To N
                        Sum = Sum + FMAS(I, J) * TMP1(J)
                      Next J
                      TMP6(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP6(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP6,ITMP)
                    For I = 1 To N
                      TMP6(I) = B(I)
                    Next I
                    For I = 1 To N
                      TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
                    Next I
                    For I = 1 To N
                      Sum = TMP3(I) + HD4 * TMP8(I)
                      For J = 1 To N
                        Sum = Sum + FMAS(I, J) * TMP1(J)
                      Next J
                      TMP7(I) = Sum
                    Next I
                    For I = 1 To N
                      B(I) = TMP7(I)
                    Next I
                    GoSub 5500 'call SOL(N,E,TMP7,ITMP)
                    For I = 1 To N
                      TMP7(I) = B(I)
                    Next I
                  End If
              End If
          End If
      Else
          If JBAND <> 0 Then
' --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX) ---
            For J = 1 To N
              'I1=IMAX(1,MUJAC+2-J)
              ia = 1: ib = MUJAC + 2 - J: GoSub 600: I1 = IMAX
              'I2=IMIN(MBJAC,N+MUJAC+1-J)
              ia = MBJAC: ib = N + MUJAC + 1 - J: GoSub 610: I2 = IMIN
              For I = I1 To I2
                E(I + MLDE, J) = -FJAC(I, J)
              Next I
              E(MDIAG, J) = E(MDIAG, J) + FAC
            Next J
            GoSub 4000 'call DECB(N,E,MLDE,MUE,ITMP,INFO)
            If INFO <> 0 Then GoTo 2080
            If IAUTNMS <> 0 Then
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
' ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
' ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS.
              For I = 1 To N
                TMP4(I) = TMP2(I)
              Next I
              For I = 1 To N
                B(I) = TMP4(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP)
              For I = 1 To N
                TMP4(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A21 * TMP4(I)
              Next I
              'call FCN(N,X,TMP1,TMP3)
              XX = X
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP5(I) = TMP3(I) + HC21 * TMP4(I)
              Next I
              For I = 1 To N
                B(I) = TMP5(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP)
              For I = 1 To N
                TMP5(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
              Next I
              'call FCN(N,X,TMP1,TMP3)
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP6(I) = TMP3(I) + HC31 * TMP4(I) + HC32 * TMP5(I)
              Next I
              For I = 1 To N
                B(I) = TMP6(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP)
              For I = 1 To N
                TMP6(I) = B(I)
              Next I
              For I = 1 To N
                TMP7(I) = TMP3(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
              Next I
              For I = 1 To N
                B(I) = TMP7(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP)
              For I = 1 To N
                TMP7(I) = B(I)
              Next I
            Else
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
' ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
' ---   3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS.
              HD1 = H * D1
              HD2 = H * D2
              HD3 = H * D3
              HD4 = H * D4
              For I = 1 To N
                TMP4(I) = TMP2(I) + HD1 * TMP8(I)
              Next I
              For I = 1 To N
                B(I) = TMP4(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP)
              For I = 1 To N
                TMP4(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A21 * TMP4(I)
              Next I
              'call FCN(N,X+C2*H,TMP1,TMP3)
              XX = X + C2 * H
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP5(I) = TMP3(I) + HD2 * TMP8(I) + HC21 * TMP4(I)
              Next I
              For I = 1 To N
                B(I) = TMP5(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP)
              For I = 1 To N
                TMP5(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
              Next I
              'call FCN(N,X+C3*H,TMP1,TMP3)
              XX = X + C3 * H
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP6(I) = TMP3(I) + HD3 * TMP8(I) + HC31 * TMP4(I) + HC32 * TMP5(I)
              Next I
              For I = 1 To N
                B(I) = TMP6(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP)
              For I = 1 To N
                TMP6(I) = B(I)
              Next I
              For I = 1 To N
                TMP7(I) = TMP3(I) + HD4 * TMP8(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
              Next I
              For I = 1 To N
                B(I) = TMP7(I)
              Next I
              GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP)
              For I = 1 To N
                TMP7(I) = B(I)
              Next I
            End If
          Else
' --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX) ---
            For J = 1 To N
              For I = 1 To N
                E(I, J) = -FJAC(I, J)
              Next I
              E(J, J) = E(J, J) + FAC
            Next J
            GoSub 4500 'call DECA(N,LDE,E,ITMP,INFO)
            If INFO <> 0 Then GoTo 2080
            If IAUTNMS <> 0 Then
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
' ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
' ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS.
              For I = 1 To N
                TMP4(I) = TMP2(I)
              Next I
              For I = 1 To N
                B(I) = TMP4(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP4,ITMP)
              For I = 1 To N
                TMP4(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A21 * TMP4(I)
              Next I
              'call FCN(N,X,TMP1,TMP3)
              XX = X
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP5(I) = TMP3(I) + HC21 * TMP4(I)
              Next I
              For I = 1 To N
                B(I) = TMP5(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP5,ITMP)
              For I = 1 To N
                TMP5(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
              Next I
              'call FCN(N,X,TMP1,TMP3)
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP6(I) = TMP3(I) + HC31 * TMP4(I) + HC32 * TMP5(I)
              Next I
              For I = 1 To N
                B(I) = TMP6(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP6,ITMP)
              For I = 1 To N
                TMP6(I) = B(I)
              Next I
              For I = 1 To N
                TMP7(I) = TMP3(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
              Next I
              For I = 1 To N
                B(I) = TMP7(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP7,ITMP)
              For I = 1 To N
                TMP7(I) = B(I)
              Next I
            Else
' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
' ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
' ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
' ---   3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS.
              HD1 = H * D1
              HD2 = H * D2
              HD3 = H * D3
              HD4 = H * D4
              For I = 1 To N
                TMP4(I) = TMP2(I) + HD1 * TMP8(I)
              Next I
              For I = 1 To N
                B(I) = TMP4(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP4,ITMP)
              For I = 1 To N
                TMP4(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A21 * TMP4(I)
              Next I
              'call FCN(N,X+C2*H,TMP1,TMP3)
              XX = X + C2 * H
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP5(I) = TMP3(I) + HD2 * TMP8(I) + HC21 * TMP4(I)
              Next I
              For I = 1 To N
                B(I) = TMP5(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP5,ITMP)
              For I = 1 To N
                TMP5(I) = B(I)
              Next I
              For I = 1 To N
                TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I)
              Next I
              'call FCN(N,X+C3*H,TMP1,TMP3)
              XX = X + C3 * H
              For I = 1 To N
                YY(I) = TMP1(I)
              Next I
              GoSub 500
              For I = 1 To N
                TMP3(I) = F(I)
              Next I
              For I = 1 To N
                TMP6(I) = TMP3(I) + HD3 * TMP8(I) + HC31 * TMP4(I) + HC32 * TMP5(I)
              Next I
              For I = 1 To N
                B(I) = TMP6(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP6,ITMP)
              For I = 1 To N
                TMP6(I) = B(I)
              Next I
              For I = 1 To N
                TMP7(I) = TMP3(I) + HD4 * TMP8(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I)
              Next I
              For I = 1 To N
                B(I) = TMP7(I)
              Next I
              GoSub 5500 'call SOL(N,E,TMP7,ITMP)
              For I = 1 To N
                TMP7(I) = B(I)
              Next I
            End If
          End If
      End If
      XNSOL = XNSOL + 4#
      XNFCN = XNFCN + 2#
' *** *** *** *** *** *** ***
'       ERROR ESTIMATION
' *** *** *** *** *** *** ***
      XNSTEP = XNSTEP + ONE
' ------------ NEW SOLUTION ---------------
      For I = 1 To N
        TMP1(I) = Y(I) + B1 * TMP4(I) + B2 * TMP5(I) + B3 * TMP6(I) + B4 * TMP7(I)
      Next I
' ------------ COMPUTE ERROR ESTIMATION ----------------
      XERR = ZERO
      For I = 1 To N
        S = E1 * TMP4(I) + E2 * TMP5(I) + E3 * TMP6(I) + E4 * TMP7(I)
        If ITOL = 0 Then
          'SK = ATOL + RTOL * XMAX(ABS(Y(I)), ABS(TMP1(I)))
          xa = Abs(Y(I)): xb = Abs(TMP1(I)): GoSub 620
          SK = ATOL + RTOL * XMAX
        Else
'         Multiple tolerances not implemented here.
'         SK = ATOL(I) + RTOL(I) * XMAX(ABS(Y(I)), ABS(TMP1(I)))
        End If
        XERR = XERR + Sqr(Abs(S / SK))
      Next I
      XERR = Sqr(XERR / N)
' --- COMPUTATION OF HNEW
' --- WE REQUIRE 0.2<=HNEW/H<=6.0
      'FAC=XMAX(FAC2,XMIN(FAC1,XERR^(0.25#/0.9#)))
      xa = FAC1: xb = XERR ^ (0.25 / 0.9): GoSub 630
      xa = FAC2: xb = XMIN: GoSub 620: FAC = XMAX
      HNEW = H / FAC
'  *** *** *** *** *** *** ***
'  IS THE XERROR SMALL ENOUGH ?
'  *** *** *** *** *** *** ***
      If XERR <= ONE Then
' --- STEP IS ACCEPTED ---
         XNACCPT = XNACCPT + ONE
         For I = 1 To N
           Y(I) = TMP1(I)
         Next I
         XXOLD = X
         X = X + H
         If IOUT <> 0 Then
           'SOLOUT(XNACCPT+1,XXOLD,X,Y,N,IRTRN)
         End If
         If IRTRN < 0 Then GoTo 2079
         If Abs(HNEW) > HMAXN Then HNEW = POSNEG * HMAXN
         If IREJECT <> 0 Then
           'HNEW=POSNEG*XMIN(ABS(HNEW),ABS(H))
           xa = Abs(HNEW): xb = Abs(H): GoSub 630
           HNEW = POSNEG * XMIN
         End If
         IREJECT = 0
         IRJECT2 = 0
         H = HNEW
         GoTo 2001
      Else
' --- STEP IS IREJECTED ---
         If IRJECT2 <> 0 Then HNEW = H * FACREJ
         If IREJECT <> 0 Then IRJECT2 = 1
         IREJECT = 1
         H = HNEW
         If XNACCPT >= ONE Then NREJCT = NREJCT + 1
         GoTo 2002
      End If
' --- EXIT SECTION ---
2080  Print " MATRIX E IS SINGULAR, INFO = "; INFO
      NSING = NSING + 1
      If NSING >= 5 Then GoTo 2079
      H = H * HALF
      GoTo 2002
2079  Print
      Print " EXIT OF ROS4 AT X="; X; "   H="; H
      IDID = -1
Return 'from RO4COR

3000  'SHAMP(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4,
            'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
         A21 = 2#
         A31 = 48# / 25#
         A32 = 6# / 25#
         C21 = -8#
         C31 = 372# / 25#
         C32 = 12# / 5#
         C41 = -112# / 125#
         C42 = -54# / 125#
         C43 = -2# / 5#
         B1 = 19# / 9#
         B2 = 1# / 2#
         B3 = 25# / 108#
         B4 = 125# / 108#
         E1 = 17# / 54#
         E2 = 7# / 36#
         E3 = 0#
         E4 = 125# / 108#
         GAMMA = HALF
         C2 = 1#
         C3 = 0.6
         D1 = 0.5
         D2 = -1.5
         D3 = 2.42
         D4 = 0.116
Return

3001  'GRK4A(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4,
            'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
        A21 = 1.10886075949367
        A31 = 2.37708526198336
        A32 = 0.185011498889969
        C21 = -4.92018840239764
        C31 = 1.05558868604858
        C32 = 3.35181726766894
        C41 = 3.84686900704931
        C42 = 3.42710924126818
        C43 = -2.16240884875326
        B1 = 1.84568324040584
        B2 = 0.13697968943605
        B3 = 0.712909778329156
        B4 = 0.632911392405063
        E1 = 4.83187017720177E-02
        E2 = -0.647110865104951
        E3 = 0.218687666050024
        E4 = -0.632911392405063
        GAMMA = 0.395
        C2 = 0.438
        C3 = 0.87
        D1 = 0.395
        D2 = -0.372672395484092
        D3 = 6.62919654457149E-02
        D4 = 0.434094696256863
Return

3002  'GRK4T(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4,
            'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
        A21 = 2#
        A31 = 4.52470820737312
        A32 = 4.16352878859765
        C21 = -5.07167533877632
        C31 = 6.02015272865079
        C32 = 0.159750684672712
        C41 = -1.85634361868611
        C42 = -8.50538085817983
        C43 = -2.08407513602319
        B1 = 3.95750374664078
        B2 = 4.62489238836331
        B3 = 0.617477263875011
        B4 = 1.28261294526904
        E1 = 2.302155402933
        E2 = 3.07363448539262
        E3 = -0.873280801804503
        E4 = -1.28261294526904
        GAMMA = 0.231
        C2 = 0.462
        C3 = 0.880208333333333
        D1 = 0.231
        D2 = -0.039629667752443
        D3 = 0.550778939578913
        D4 = -5.53509845705276E-02
Return

3003  'VELDS(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4,
            'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
' --- METHOD GIVEN BY VAN VELDHUIZEN ---
        A21 = 2#
        A31 = 1.75
        A32 = 0.25
        C21 = -8#
        C31 = -8#
        C32 = -1#
        C41 = 0.5
        C42 = -0.5
        C43 = 2#
        B1 = 1.33333333333333
        B2 = 0.666666666666667
        B3 = -1.33333333333333
        B4 = 1.33333333333333
        E1 = -0.333333333333333
        E2 = -0.333333333333333
        E3 = -0#
        E4 = -1.33333333333333
        GAMMA = 0.5
        C2 = 1#
        C3 = 0.5
        D1 = 0.5
        D2 = -1.5
        D3 = -0.75
        D4 = 0.25
Return

3004 'VELDD(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4,
           'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
' --- METHOD GIVEN BY VAN VELDHUIZEN ---
        A21 = 2#
        A31 = 4.81223436269544
        A32 = 4.57814695674784
        C21 = -5.33333333333333
        C31 = 6.10052967884825
        C32 = 1.80473679737843
        C41 = -2.54051545663475
        C42 = -9.44374632891521
        C43 = -1.98847175321599
        B1 = 4.28933925465454
        B2 = 5.03609848285141
        B3 = 0.608573642067392
        B4 = 1.35595894120115
        E1 = 2.17567278753176
        E2 = 2.95091122257574
        E3 = -0.785974454488743
        E4 = -1.35595894120115
        GAMMA = 0.225708114822568
        C2 = 0.451416229645136
        C3 = 0.875592894601846
        D1 = 0.225708114822568
        D2 = -4.59940350268058E-02
        D3 = 0.517759050494408
        D4 = -3.80562393805443E-02
Return

3005  'LSTAB(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4,
            'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
' --- AN L-STABLDE METHOD ---
        A21 = 2#
        A31 = 1.86794363780392
        A32 = 0.234444971139916
        C21 = -7.13761503641231
        C31 = 2.58070808795146
        C32 = 0.651595007644798
        C41 = -2.13714899438253
        C42 = -0.321466969123763
        C43 = -0.694974250178178
        B1 = 2.25557007341874
        B2 = 0.287049326218679
        B3 = 0.435317943184018
        B4 = 1.09350225240916
        E1 = -0.281543193214115
        E2 = -7.27619912493892E-02
        E3 = -0.108219620149531
        E4 = -1.09350225240916
        GAMMA = 0.57282
        C2 = 1.14564
        C3 = 0.65521686381559
        D1 = 0.57282
        D2 = -1.76919389131923
        D3 = 0.759263343792048
        D4 = -0.104902108710045
Return

4000 'Subroutine DECB(N, E, MLDE, MUE, ITMP, INFO)
' Labels 4007, 4020, 4030, 4045, 4060, 4070, 4080
'-----------------------------------------------------------------------
'  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
'  MATRIX WITH LOWER BANDWIDTH MLDE AND UPPER BANDWIDTH MUE
'  INPUTS:
'     N       ORDER OF THE ORIGINAL MATRIX A.
'     E       CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
'             OF THE MATRIX ARE STORED IN THE COLUMNS OF  E  AND
'             THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
'             MLDE+1 THROUGH 2*MLDE+MUE+1 OF E.
'     MLDE    LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
'     MUE     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
'  OUTPUTS:
'     E       AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
'             THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
'     ITMP    INDEX VECTOR OF PIVOT INDICES.
'     ITMP(N) (-1)^(NUMBER OF INTERCHANGES) OR O.
'     INFO     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
'              SINGULAR AT STAGE K.
'  NOTE: USE SOLB TO OBTAIN SOLUTION OF LINEAR SYSTEM.
'  DETERM(E) = ITMP(N)*E(MD,1)*E(MD,2)*...*E(MD,N)  WITH MD=MLDE+MUE+1.
'  IF ITMP(N)=O, E IS SINGULAR, SOLB WILL DIVIDE BY ZERO.
'
'  REFERENCE..
'     THIS IS A MODIFICATION OF
'     C. B. MOLDER, ALGORITHM 423, LINEAR EQUATION SOLVER,
'     C.A.C.M. 15 (1972), P. 274.
'-----------------------------------------------------------------------
      INFO = 0
      ITMP(N) = 1
      MD = MLDE + MUE + 1
      MD1 = MD + 1
      JU = 0
      If MLDE = 0 Then GoTo 4070
      If N = 1 Then GoTo 4070
      If N < MUE + 2 Then GoTo 4007
      For J = MUE + 2 To N
        For I = 1 To MLDE
          E(I, J) = ZERO
        Next I
      Next J
4007  NM1 = N - 1
      For K = 1 To NM1
        KP1 = K + 1
        M = MD
        'MDL = IMIN(MLDE,N-K) + MD
        ia = MLDE: ib = N - K: GoSub 610
        MDL = IMIN + MD
        For I = MD1 To MDL
          If Abs(E(I, K)) > Abs(E(M, K)) Then M = I
        Next I
        ITMP(K) = M + K - MD
        t = E(M, K)
        If M = MD Then GoTo 4020
        ITMP(N) = -ITMP(N)
        E(M, K) = E(MD, K)
        E(MD, K) = t
4020    If t = ZERO Then GoTo 4080
        t = ONE / t
        For I = MD1 To MDL
          E(I, K) = -E(I, K) * t
        Next I
        'JU = IMIN(IMAX(JU,MUE+ITMP(K)),N)
        ia = JU: ib = MUE + ITMP(K): GoSub 600
        ia = IMAX: ib = N: GoSub 610: JU = IMIN
        MM = MD
        If JU < KP1 Then GoTo 4060
        For J = KP1 To JU
          M = M - 1
          MM = MM - 1
          t = E(M, J)
          If M = MM Then GoTo 4030
          E(M, J) = E(MM, J)
          E(MM, J) = t
4030      If t = ZERO Then GoTo 4045
          JK = J - K
          For I = MD1 To MDL
            IJK = I - JK
            E(IJK, J) = E(IJK, J) + E(I, K) * t
          Next I
4045    Next J
4060  Next K
4070  K = N
      If E(MD, N) = ZERO Then GoTo 4080
      Return
4080  INFO = K
      ITMP(N) = 0
Return 'DECB

4500 'Subroutine DECA(N, E, ITMP, INFO)
' REAL DOUBLE PRECISION VERSION OF DEC
' Labels: 4520, 4550, 4570, 4580
'-----------------------------------------------------------------------
'  GENERAL MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
'  INPUT..
'     N = ORDER OF MATRIX.
'     E = MATRIX TO BE TRIANGULARIZED.
'  OUTPUT..
'     E(I,J), I <= J = UPPER TRIANGULAR FACTOR, U.
'     E(I,J), I > J  = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
'     ITMP(K), K < N = INDEX OF K-TH PIVOT ROW.
'     ITMP(N) = (-1)^(NUMBER OF INTERCHANGES) OR O.
'     INFO = 0 IF MATRIX E IS NONSINGULAR, OR K IF FOUND TO BE
'            SINGULAR AT STAGE K.
'  NOTE: USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM.
'  DETERM(E) = ITMP(N)*E(1,1)*E(2,2)*...*E(N,N).
'  IF ITMP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
'
'  REFERENCE..
'     C. B. MOLDER, ALGORITHM 423, LINEAR EQUATION SOLVER,
'     C.A.C.M. 15 (1972), P. 274.
'------------------------------------------------------------------------
      INFO = 0
      ITMP(N) = 1
      If N = 1 Then GoTo 4570
      NM1 = N - 1
      For K = 1 To NM1
        KP1 = K + 1
        M = K
        For I = KP1 To N
          If Abs(E(I, K)) > Abs(E(M, K)) Then M = I
        Next I
        ITMP(K) = M
        t = E(M, K)
        If M = K Then GoTo 4520
        ITMP(N) = -ITMP(N)
        E(M, K) = E(K, K)
        E(K, K) = t
4520    If t = ZERO Then GoTo 4580
        t = ONE / t
        For I = KP1 To N
          E(I, K) = -E(I, K) * t
        Next I
        For J = KP1 To N
          t = E(M, J)
          E(M, J) = E(K, J)
          E(K, J) = t
          If t = ZERO Then GoTo 4550
          For I = KP1 To N
            E(I, J) = E(I, J) + E(I, K) * t
          Next I
4550    Next J
      Next K
4570  K = N
      If E(N, N) = ZERO Then GoTo 4580
      Return
4580  INFO = K
      ITMP(N) = 0
Return


5000 'Subroutine SOLB(N, E, MLDE, MUE, B, ITMP)
' Labels 5025, 5050
'-----------------------------------------------------------------------
'  SOLUTION OF A BANDED LINEAR SYSTEM, E*X = B.
'  INPUTS:
'    N      ORDER OF MATRIX A.
'    E      TRIANGULARIZED MATRIX OBTAINED FROM DECB.
'    MLDE   LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
'    MUE    UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
'    B      RIGHT HAND SIDE VECTOR.
'    ITMP   PIVOT VECTOR OBTAINED FROM DECB.
'  DO NOT USE IF DECB HAS SET INFO <> 0.
'  OUTPUT:
'    B      SOLUTION VECTOR.
'-----------------------------------------------------------------------
      MD = MLDE + MUE + 1
      MD1 = MD + 1
      MDM = MD - 1
      NM1 = N - 1
      If MLDE = 0 Then GoTo 5025
      If N = 1 Then GoTo 5050
      For K = 1 To NM1
        M = ITMP(K)
        t = B(M)
        B(M) = B(K)
        B(K) = t
        'MDL = IMIN(MLDE,N-K) + MD
        ia = MLDE: ib = N - K: GoSub 610
        MDL = IMIN + MD
        For I = MD1 To MDL
          IMD = I + K - MD
          B(IMD) = B(IMD) + E(I, K) * t
        Next I
      Next K
5025  For KB = 1 To NM1
        K = N + 1 - KB
        B(K) = B(K) / E(MD, K)
        t = -B(K)
        KMD = MD - K
        'LM = IMAX(1,KMD+1)
        ia = 1: ib = KMD + 1: GoSub 600
        LM = IMAX
        For I = LM To MDM
          IMD = I - KMD
          B(IMD) = B(IMD) + E(I, K) * t
        Next I
      Next KB
5050  B(1) = B(1) / E(MD, 1)
Return

5500  'Subroutine SOL (N, E, B, ITMP)
' Label: 5550
'-----------------------------------------------------------------------
'  SOLUTION OF A GENERAL LINEAR SYSTEM, E*X = B.
'  INPUTS:
'    N = ORDER OF MATRIX.
'    E = TRIANGULARIZED MATRIX OBTAINED FROM DECA.
'    B = RIGHT HAND SIDE VECTOR.
'    ITMP = PIVOT VECTOR OBTAINED FROM DEC.
'  DO NOT USE IF DEC HAS SET INFO <> 0.
'  OUTPUT:
'    B = SOLUTION VECTOR.
'-----------------------------------------------------------------------
      If N = 1 Then GoTo 5550
      NM1 = N - 1
      For K = 1 To NM1
        KP1 = K + 1
        M = ITMP(K)
        t = B(M)
        B(M) = B(K)
        B(K) = t
        For I = KP1 To N
          B(I) = B(I) + E(I, K) * t
        Next I
      Next K
      For KB = 1 To NM1
        KM1 = N - KB
        K = KM1 + 1
        B(K) = B(K) / E(K, K)
        t = -B(K)
        For I = 1 To KM1
          B(I) = B(I) + E(I, K) * t
        Next I
      Next KB
5550  B(1) = B(1) / E(1, 1)
Return

'end of file tros4.bas