1SECTION : 5.2.1.1.1.1        (DECEMBER 1979)                     PAGE 1
 
 
 
 SECTION 5.2.1.1.1.1 CONTAINS NINE PROCEDURES FOR INTIAL-VALUE PROBLEMS
 FOR FIRST ORDER  ORDINARY DIFERENTIAL EQUATIONS.
 
 A.  RK1 SOLVES AN IVP FOR A SINGLE ODE BY MEANS OF A
     5-TH ORDER RUNGE-KUTTA METHOD.
 B.  RKE SOLVES AN IVP FOR A SYSTEM OF ODE'S BY MEANS OF A
     5-TH ORDER RUNGE-KUTTA METHOD.
 
 C.  RK4A SOLVES AN IVP FOR A SINGLE ODE BY MEANS OF A 5-TH ORDER
     RUNGE-KUTTA METHOD. RK4A INTERCHANGES THE DEPENDENT AND
     INDEPENDENT VARIABLE.
 
 D.  RK4NA SOLVES AN IVP FOR A SYSTEM OF ODE'S BY MEANS OF A
     5-TH ORDER RUNGE-KUTTA METHOD. RK4NA INTERCHANGES THE
     DEPENDENT AND ONE INDEPENDENT VARIABLE.
 
 E.  RK5NA SOLVES AN IVP FOR A SYSTEM OF ODE'S BY MEANS OF A
     5-TH ORDER RUNGE-KUTTA METHOD. RK5NA USES THE ARC LENGTH
     AS INTEGRATION VARIABLE.
 
 F.  MULTISTEP  SOLVES  AN IVP FOR A SYSTEM OF ODE'S BY MEANS OF A
     LINEAR MULTISTEP METHOD. IT USES EITHER THE BACKWARD
     DIFFERENTIATION METHODS, OR THE ADAMS-BASHFORTH-MOULTON-METHOD.
 
 G.  DIFFSYS  SOLVES AN IVP FOR A SYSTEM OF ODE'S BY MEANS OF
     A HIGH ORDER EXTRAPOLATION METHOD BASED ON THE MODIFIED
     MIDPOINT RULE.
 
 H.  ARK SOLVES AN IVP FOR A LARGE SYSTEM OF ODE'S WHICH IS OBTAINED
     FROM SEMI-DISCRETIZATION OF AN INITIAL BOUNDARY VALUE PROBLEM
     FOR A PARABOLIC OR HYPERBOLIC EQUATION. ARK IS BASED ON
     STABILIZED, EXPICIT RUNGE-KUTTA METHODS OF LOW ORDER
 
 I.  EFRK  SOLVES  AN IVP FOR A SYSTEM OF ODE'S BY MEANS OF AN
     EXPONENTIALLY FITTED EXPLICIT RUNGE-KUTTA METHOD OF FIRST,
     SECOND OR THIRD ORDER.
 
 
 
 RK1 AND RKE ARE INTENDED FOR NON-STIFF EQUATIONS,WHILE RK4A, RK4NA
 AND RK5NA   ARE INTENDED FOR NON-SIFF EQUATIONS WHERE DERIVATIVES
 BECOME VERY LARGE, SUCH AS IN THE NEIGHBOURHOOD OF SINGULARITIES.
 MULTISTEP CAN BE USED FOR NON-STIF, AS WELL AS STIFF EQUATIONS.
 DIFFSYS SHOULD BE USED FOR PROBLEMS FOOR WHICH A VERY HIGH ACCURACY
 IS DESIRED. ARK IS RECOMMENDED FOR THE INTEGRATION OF SEMI-DISCRETE
 PARABOLIC OR HYPERBOLIC PROBLEMS. EFRK IS A SPECIAL PURPOSE PROCEDURE
 FOR STIFF EQUATIONS WITH A KNOWN, CLUSTERED EIGENVALUE SPECTRUM.
 EXCEPT EFRK, ALL PROCEDURES PERFORM STEPSIZE CONTROL.
1SECTION : 5.2.1.1.1.1.A      (FEBRUARY 1979)                     PAGE 1
 
 
 
 PROCEDURE : RK1
 
 
 AUTHOR:J.A.ZONNEVELD.
 
 
 CONTRIBUTORS: M.BAKKER AND I.BRINK.
 
 
 INSTITUTE: MATHEMATICAL CENTRE.
 
 
 RECEIVED: 730715.
 
 
 BRIEF DESCRIPTION:
 
     RK1 SOLVES AN INITIAL VALUE PROBLEM FOR A SINGLE FIRST ORDER
     ORDINARY DIFFERENTIAL EQUATION    DY / DX  = F(X,Y).
     THE EQUATION IS SUPOSED TO BE NON-STIFF.
 
 
 KEYWORDS:
 
     INITIAL VALUE PROBLEM.
     SINGLE FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK1(X,A,B,Y,YA,FXY,E,D,FI);
     "VALUE" B,FI;
     "REAL" X,A,B,Y,YA,FXY;
     "BOOLEAN" FI;
     "ARRAY" E,D;
     "CODE" 33010;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     X:    <VARIABLE>;
           THE INDEPENDENT VARIABLE.
           UPON COMPLETION OF A CALL
           X IS EQUAL TO B;
     A:    <ARITHMETIC EXPRESSION>;
           ENTRY: THE INITIAL VALUE OF X;
     B:    <ARITHMETIC EXPRESSION>;
           ENTRY: A VALUE PARAMETER,GIVING THE END VALUE OF X;
     Y:    <VARIABLE>;
           THE DEPENDENT VARIABLE;
1SECTION : 5.2.1.1.1.1.A      (FEBRUARY 1979)                     PAGE 2
 
 
 
     YA:   <ARITHMETIC EXPRESSION>;
           ENTRY :  THE VALUE OF Y AT X=A;
     FXY:  <ARITHMETIC EXPRESSION>;
           AN EXPRESSION,DEPENDING ON X AND Y,GIVING THE VALUE OF DY/DX;
     E:    <ARRAY IDENTIFIER>;
           "ARRAY" E[1:2];
           ENTRY:
           E[1] : A RELATIVE TOLERANCE,
           E[2] : AN ABSOLUTE TOLERANCE ;
     D:    <ARRAY IDENTIFIER>;
           "ARRAY" D[1:4];
           EXIT:
           ENTIER(D[1]+.5) GIVES THE NUMBER OF STEPS SKIPPED;
           D[2] : EQUALS THE STEP LENGTH;
           D[3] : IS EQUAL TO B;
           D[4] : IS EQUAL TO Y(B);
     FI:   <BOOLEAN EXPRESSION>;
           IF FI="TRUE" RK1 INTEGRATES FROM X=A TO X=B WITH INITIAL VALUE
           VALUE Y(A)=YA AND TRIAL STEP B-A. IF  FI="FALSE" RK1 INTEGRATES.
           FROM X=D[3] TO X=B WITH INITIAL VALUE Y(D[3])=D[4] ANSD FIRST
           1TEP H=D[2]*SIGN(B-D[3]), WHILE A AND YA ARE IGNORED.
 
 
 PROCEDURES USED: NONE.
 
 
 REQUIRED CENTRAL MEMORY : NEGLIGIBLE SMALL.
 
 
 METHOD AND PERFORMANCE:
 RK1 IS BASED ON AN EXPLICIT, 5-TH ORDER RUNGE-KUTTA METHOD AND
 IS PROVIDED WITH STEPLENGTH AND ERRORCONTROL. THE ERRORCONTROL
 IS BASED ON THE LAST TERM OF THE TAYLOR SERIES WICH IS TAKEN INTO
 ACCOUNT. A STEP IS REJECTED IF THE ABSOLUTE VALUE OF THIS LAST TERM
 GREATER THEN  (ABS(FXY)*E[1]+E[2])*ABS(H)/INT, WHERE
 INT = ABS( B-( "IF" FI "THEN" A "ELSE" D[3])) DENOTES THE LENGTH
 OF THE INTEGRATION INTERVAL, OTHERWISE, A STEP IS ACCEPTED.
 RK1 USES AS ITS MINIMAL ABBSOLUTE STEPLENGTH  HMIN = E[1]*INT+E[2].
 IF A STEP OF LENGTH ABS(H) = HMIN  IS REJECTED,  THE STEP IS SKIPPED.
 FOR FURTHER DETAILS SEE [1].
1SECTION : 5.2.1.1.1.1.A      (FEBRUARY 1979)                     PAGE 3
 
 
 
 
 REFERENCES:
     [1]J.A.ZONNEVELD.
          AUTOMATIC NUMERICAL INTEGRATION.
          MATHEMATICAL CENTRE TRACT 8(1970).
 
 
 EXAMPLE OF USE:
 
     THE SOLUTION OF THE DIFFERENTIAL EQUATION DY/DX=-Y
     WITH INITIAL CONDITION Y(0)=1 AT X=1 IS COMPUTED
     BY MEANS OF THE FOLLOWING PROGRAM:
 
     "BEGIN"
     "REAL" X,Y;
     "BOOLEAN" FI,FIRST;
     "REAL" "ARRAY" E[1:2],D[1:4];
 
     "PROCEDURE" RK1(X,A,B,Y,YA,FXY,E,D,FI); "CODE" 33010;
 
     E[1]:=+"-4;E[2]:=+"-4;FIRST:="TRUE";
     RK1(X,0,1,Y,1,-Y,E,D,FIRST);
     OUTPUT(61,"("//10B"("X=")".12D"2D,//10B"("Y=")".12D"2D,
     10B"("YEXACT=")".12D"2D")",X,Y,EXP(-X));
     "END"
         "EOP"
 
     IT DELIVERS WITH E[1]=E[2]="-4:
           X=.100000000000"01
 
 
           Y=.367876846355"00          YEXACT=.367879441171"00
 
 SOURCE TEXT(S):
0"CODE"" 33010;
    "PROCEDURE" RK1(X, A, B, Y, YA, FXY, E, D, FI);
    "VALUE" B, FI; "REAL" X, A, B, Y, YA, FXY; "BOOLEAN" FI;
    "ARRAY" E, D;
    "BEGIN" "REAL" E1, E2, XL, YL, H, INT, HMIN, ABSH, K0, K1,
       K2, K3, K4, K5, DISCR, TOL, MU, MU1, FH, HL;
       "BOOLEAN" LAST, FIRST, REJECT;                          "COMMENT"
1SECTION : 5.2.1.1.1.1.A      (AUGUST 1974)                       PAGE 4
                                                                  ;
 
 
       "IF" FI "THEN"
       "BEGIN" D[3]:= A; D[4]:= YA "END";
       D[1]:= 0; XL:= D[3]; YL:= D[4];
       "IF" FI "THEN" D[2]:= B - D[3]; ABSH:= H:= ABS(D[2]);
       "IF" B - XL < 0 "THEN" H:= - H; INT:= ABS(B - XL);
       HMIN:= INT * E[1] + E[2]; E1:= E[1] / INT;
       E2:= E[2] / INT; FIRST:= "TRUE"; "IF" FI "THEN"
       "BEGIN" LAST:= "TRUE"; "GOTO" STEP "END";
    TEST: ABSH:= ABS(H); "IF" ABSH < HMIN "THEN"
       "BEGIN" H:= "IF" H > 0 "THEN" HMIN "ELSE" - HMIN; ABSH:= HMIN
       "END";
       "IF" H >= B - XL "EQUIV" H >= 0 "THEN"
       "BEGIN" D[2]:= H; LAST:= "TRUE"; H:= B - XL;
          ABSH:= ABS(H)
       "END"
       "ELSE" LAST:= "FALSE";
    STEP: X:= XL; Y:= YL; K0:= FXY * H;
       X:= XL + H / 4.5; Y:= YL + K0 / 4.5;
       K1:= FXY * H; X:= XL + H / 3;
       Y:= YL + (K0 + K1 * 3) / 12; K2:= FXY * H;
       X:= XL + H * .5; Y:= YL + (K0 + K2 * 3) / 8;
       K3:= FXY * H; X:= XL + H * .8;
       Y:= YL + (K0 * 53 - K1 * 135 + K2 * 126 + K3 * 56)
       / 125; K4:= FXY * H; X:= "IF" LAST "THEN" B "ELSE" XL + H;
       Y:= YL + (K0 * 133 - K1 * 378 + K2 * 276 + K3 * 112
       + K4 * 25) / 168; K5:= FXY * H;
       DISCR:= ABS(K0 * 21 - K2 * 162 + K3 * 224 - K4 * 125
       + K5 * 42) / 14; TOL:= ABS(K0) * E1 + ABSH * E2;
       REJECT:= DISCR > TOL; MU:= TOL / (TOL + DISCR) + .45;
       "IF" REJECT "THEN"
       "BEGIN" "IF" ABSH <= HMIN "THEN"
          "BEGIN" D[1]:= D[1] + 1; Y:= YL; FIRST:= "TRUE";
             "GOTO" NEXT
          "END";
          H:= MU * H; "GOTO" TEST
       "END";
       "IF" FIRST "THEN"
       "BEGIN" FIRST:= "FALSE"; HL:= H; H:= MU * H; "GOTO" ACC
       "END";
       FH:= MU * H / HL + MU - MU1; HL:= H; H:= FH * H;
    ACC: MU1:= MU;
       Y:= YL + ( - K0 * 63 + K1 * 189 - K2 * 36 - K3 * 112
       + K4 * 50) / 28; K5:= FXY * HL;
       Y:= YL + (K0 * 35 + K2 * 162 + K4 * 125 + K5 * 14)
       / 336;
 
    NEXT: "IF" B ^= X "THEN"
       "BEGIN" XL:= X; YL:= Y; "GOTO" TEST "END";
       "IF" "NOT"LAST "THEN" D[2]:= H; D[3]:= X; D[4]:= Y
    "END" RK1;
         "EOP"
1SECTION : 5.2.1.1.1.1.B     (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : RKE. 
        
        
 AUTHOR: P.A. BEENTJES.     
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 740520.
        
        
 BRIEF DESCRIPTION: 
        
     RKE  SOLVES AN INITIAL VALUE PROBLEM FOR A SYSTEM OF FIRST     
     ORDER ORDINARY DIFFERENTIAL EQUATIONS  DY / DX = F(X,Y).       
     THE SYSTEM IS SUPPOSED TO BE NON-STIFF.    
        
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEM. 
     SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS. 
        
1SECTION : 5.2.1.1.1.1.B     (FEBRUARY 1979)                     PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
        
     "PROCEDURE" RKE (X, XE, N, Y, DER, DATA, FI, OUT);   
     "VALUE" N, FI; "INTEGER" N; "REAL" X, XE; "BOOLEAN" FI;
     "ARRAY" Y, DATA;       
     "PROCEDURE" DER, OUT;  
     "CODE" 33033;
        
     RKE  INTEGRATES  THE  SYSTEM OF ORDINARY DIFFERENTIAL  EQUATIONS 
      DY / DX = F(X, Y),  FROM X = X0 TO X = XE WHERE Y(X0) = Y0.   
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:      <VARIABLE>;    
             THE INDEPENDENT VARIABLE;
             ENTRY: THE INITIAL VALUE X0;       
     XE:     <ARITHMETIC EXPRESSION>; 
             THE FINAL VALUE OF X;    
     N:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF EQUATIONS OF THE SYSTEM;       
     Y:      <ARRAY IDENTIFIER>;      
             "ARRAY" Y[1 : N];
             THE DEPENDENT VARIABLES; 
             ENTRY: THE INITIAL VALUES AT X = X0; 
             EXIT : THE VALUES OF THE SOLUTION AT X = XE; 
     DER:    <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS: 
             "PROCEDURE" DER (T, V); "VALUE" T; "REAL" T; "ARRAY" V;
             THIS PROCEDURE  PERFORMS  AN  EVALUATION OF THE RIGHT HAND       
             SIDE  OF THE SYSTEM WITH DEPENDENT VARIABLES V[1 : N]  AND       
             INDEPENDENT  VARIABLE  T; UPON  COMPLETION  OF DER     
             THE  RIGHT HAND SIDE SHOULD BE  OVERWRITTEN  ON  V[1 : N];       
     DATA:   <ARRAY IDENTIFIER>;      
             "ARRAY" DATA[1 : 6];     
             IN ARRAY DATA ONE SHOULD GIVE:     
                 DATA[1]: THE RELATIVE TOLERANCE; 
                 DATA[2]: THE ABSOLUTE TOLERANCE; 
             AFTER  EACH  STEP DATA[3:6] CONTAINS :       
                 DATA[3]: THE   STEPLENGTH   USED  FOR THE  LAST  STEP;       
                 DATA[4]: THE  NUMBER  OF INTEGRATION STEPS  PERFORMED;       
                 DATA[5]: THE  NUMBER  OF  INTEGRATION STEPS  REJECTED;       
                 DATA[6]: THE  NUMBER  OF  INTEGRATION  STEPS  SKIPPED;       
                          IF  UPON  COMPLETION  OF   RKE  DATA[6] > 0 ,       
                          RESULTS SHOULD BE CONSIDERED MOST CRITICALLY;       
     FI:     <BOOLEAN EXPRESSION>;    
             IF  FI = "TRUE"  THE  INTEGRATION  STARTS  AT  X0  WITH  A       
             TRIAL  STEP  XE - X0; IF  FI = "FALSE"  THE INTEGRATION IS       
             CONTINUED   WITH   A  STEPLENGTH  DATA[3] * SIGN(XE - X0);       
     OUT:    <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS: 
             "PROCEDURE" OUT; 
             AFTER  EACH INTEGRATION STEP PERFORMED OUT CAN BE USED TO
             OBTAIN INFORMATION FROM THE SOLUTION PROCES, E.G. THE VALUES     
             OF X, Y[1 : N] AND  DATA[3 : 6]. POUT CAN ALSO BE USED TO
             UPDATE DATA.   
        
1SECTION : 5.2.1.1.1.1.B     (FEBRUARY 1979)                     PAGE 3      
        
        
        
 PROCEDURES USED : NONE.    
        
        
 REQUIRED CENTRAL MEMORY:   
        
     CIRCA  5 * N MEMORY PLACES.      
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD UPON WHICH THE PROCEDUDRE IS BASED, IS A MEMBER OF A
     CLASS OF FIFTH ORDER  RUNGE KUTTA  FORMULAS  PRESENTED IN REFERENCE      
     [1]. AUTOMATIC STEPSIZE CONTROL IS IMPLEMENTED IN A WAY AS  PROPOSED     
     IN REFERENCE [2]. FOR TESTRESULTS AND FURTHER INFORMATION      
     SEE  REFERENCE [3].    
        
        
 REFERENCES:      
        
     [1]. R. ENGLAND.       
             ERROR ESTIMATES FOR RUNGE KUTTA TYPE SOLUTIONS TO  SYSTEMS       
             OF ORDINARY DIFFERENTIAL EQUATIONS.
             THE COMPUTER JOURNAL , VOLUME 12, P 166 - 169, 1969.   
        
     [2]. J.A. ZONNEVELD.   
             AUTOMATIC NUMERICAL INTEGRATION.   
             MATH. CENTRE TRACT 8(1970).
        
     [3]. P.A. BEENTJES.    
             SOME  SPECIAL FORMULAS OF THE ENGLAND CLASS OF FIFTH ORDER       
             RUNGE - KUTTA SCHEMES.   
             MATH. CENTRE REPORT NW 24/74.      
        
1SECTION : 5.2.1.1.1.1.B     (DECEMBER 1975)                     PAGE 4      
        
        
        
 EXAMPLE OF USE:  
        
     THE SOLUTION AT T = 1 AND T = -1 OF THE SYSTEM       
        
     DX / DT = Y - Z,       
     DY / DT = X * X + 2 * Y + 4 * T, 
     DZ / DT = X * X + 5 * X + Z * Z + 4 * T,   
        
     WITH  X = Y = 0  AND  Z = 2  AT  T = 0,    
        
     CAN BE OBTAINED BY THE FOLLOWING PROGRAM:  
        
        
     "BEGIN" "REAL" T, TE; "ARRAY" Y[1 : 3], DATA[1 : 6]; 
        
        "PROCEDURE" RKE(X, XE, N, Y, DER, DATA, FI, OUT); 
        "CODE" 33033;       
        
        "PROCEDURE" RHS(T, Y); "VALUE" T; "REAL" T; "ARRAY" Y;      
        "BEGIN" "REAL" XX, YY, ZZ;    
           XX:= Y[1]; YY:= Y[2]; ZZ:= Y[3];     
           Y[1]:= YY - ZZ;  
           Y[2]:= XX * XX + 2 * YY + 4 * T;     
           Y[3]:= XX * (XX + 5) + 2 * ZZ + 4 * T
        "END" RHS;
        
        "PROCEDURE" INFO;   
        "IF" T = TE "THEN"  
        "BEGIN" "REAL" ET, T2, AEX, AEY, AEZ, REX, REY, REZ;
           ET:= EXP(T); T2:= 2 * T;   
           REX:= -ET * SIN(T2); AEX:= REX - Y[1]; REX:= ABS(AEX / REX);       
           REY:= ET * ET * (8 + 2 * T2 - SIN(2 * T2)) / 8 - T2 - 1; 
           REZ:= ET * (SIN(T2) + 2 * COS(T2)) + REY;      
           AEY:= REY - Y[2]; REY:= ABS(AEY / REY); AEZ:= REZ - Y[3];
           REZ:= ABS(AEZ / REZ);      
           OUTPUT(61, "(""(" T = ")", +D, //,   
           "(" RELATIVE AND ABSOLUTE ERRORS IN X, Y AND Z :")", //, 
           "("   RE(X)   RE(Y)   RE(Z)   AE(X)   AE(Y)   AE(Z)")", //,
           6(B,-.2D"+D), //,
           "(" NUMBER OF INTEGRATION STEPS PERFORMED :")",4ZD,/,    
           "(" NUMBER OF INTEGRATION STEPS SKIPPED   :")",4ZD,/,    
           "(" NUMBER OF INTEGRATION STEPS REJECTED  :")",4ZD,///")", 
           T, REX, REY, REZ, ABS(AEX), ABS(AEY), ABS(AEZ),
           DATA[4], DATA[6], DATA[5]) 
        "END" INFO; 
        
        TE:= 1;   
     LEFT:
        Y[1]:= Y[2]:= 0; Y[3]:= 2; T:=0;
        DATA[1]:= DATA[2]:= "-5;      
        RKE(T, TE, 3, Y, RHS, DATA, "TRUE", INFO);
        "IF" TE = 1 "THEN" "BEGIN" TE:= -1; "GOTO" LEFT "END"       
     "END"
        
        
1SECTION : 5.2.1.1.1.1.B     (DECEMBER 1975)                     PAGE 5      
        
        
        
     THIS PROGRAM DELIVERS: 
        
        
     T = +1       
        
     RELATIVE AND ABSOLUTE ERRORS IN X, Y AND Z : 
        
       RE(X)   RE(Y)   RE(Z)   AE(X)   AE(Y)   AE(Z)      
        
      .37"-6  .15"-5  .13"-5  .91"-6  .13"-4  .11"-4      
        
     NUMBER OF INTEGRATION STEPS PERFORMED :    9 
     NUMBER OF INTEGRATION STEPS SKIPPED   :    0 
     NUMBER OF INTEGRATION STEPS REJECTED  :    5 
        
        
     T = -1       
        
     RELATIVE AND ABSOLUTE ERRORS IN X, Y AND Z : 
        
       RE(X)   RE(Y)   RE(Z)   AE(X)   AE(Y)   AE(Z)      
        
      .22"-6  .52"-7  .19"-6  .75"-7  .55"-7  .77"-7      
        
     NUMBER OF INTEGRATION STEPS PERFORMED :   10 
     NUMBER OF INTEGRATION STEPS SKIPPED   :    0 
     NUMBER OF INTEGRATION STEPS REJECTED  :    7 
        
        
 SOURCE TEXT(S):  
0"CODE" 33033;    
    "PROCEDURE" RKE (X, XE, N, Y, DER, DATA, FI, OUT);    
    "VALUE"  FI, N; "INTEGER" N; "REAL" X, XE;  
    "BOOLEAN" FI; "ARRAY" Y, DATA;    
    "PROCEDURE" DER, OUT;   
    "BEGIN" "INTEGER" J;    
       "REAL" XT, H, HMIN, INT, HL, HT, ABSH, FHM, DISCR, TOL, MU,  
       MU1, FH, E1, E2;     
       "BOOLEAN" LAST, FIRST, REJECT; 
       "ARRAY" K0, K1, K2, K3, K4[1:N]; 
       "IF" FI "THEN"       
       "BEGIN" DATA[3]:= XE - X; DATA[4]:= DATA[5]:= DATA[6]:= 0 "END";       
       ABSH:= H:= ABS(DATA[3]);       
       "IF" XE < X "THEN" H:= - H; INT:= ABS(XE - X);     
       HMIN:= INT * DATA[1] + DATA[2];
       E1:= 12 * DATA[1] / INT; E2:= 12 * DATA[2] / INT;  
       FIRST:= "TRUE"; REJECT:= "FALSE"; "IF" FI "THEN"   
       "BEGIN" LAST:= "TRUE"; "GOTO" STEP "END";
    TEST: ABSH:= ABS(H); "IF" ABSH < HMIN "THEN"
       "BEGIN" H:= SIGN (XE - X) * HMIN; ABSH:= HMIN "END"; 
       "IF" H >= XE - X "EQUIV" H >= 0 "THEN"   
       "BEGIN" LAST:= "TRUE"; H:= XE - X; ABSH:= ABS(H) "END"       
       "ELSE" LAST:= "FALSE"; 
                                                               "COMMENT"      
1SECTION : 5.2.1.1.1.1.B     (DECEMBER 1975)                     PAGE 6      
                                                                  ; 
        
        
    STEP: "IF" ^ REJECT "THEN"
       "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K0[J]:= Y[J];    
           DER(X, K0)       
       "END";     
       HT:= .184262134833347 * H; XT:= X + HT;  
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K1[J]:= K0[J] * HT + Y[J]; 
       DER(XT, K1); 
       HT:= .690983005625053"-1 * H; XT:= 4 * HT + X;     
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K2[J]:=
       (3 * K1[J] + K0[J]) * HT + Y[J]; 
       DER(XT, K2); 
       XT:= .5 * H + X; HT:= .1875 * H; 
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K3[J]:=((1.74535599249993
       * K2[J] - K1[J]) * 2.23606797749979 + K0[J]) * HT + Y[J];    
       DER(XT, K3); 
       XT:= .723606797749979 * H + X; HT:= .4 * H;
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K4[J]:= (((.517595468166681
       * K0[J] - K1[J]) * .927050983124840 + K2[J]) * 1.46352549156242
       + K3[J]) * HT + Y[J];
       DER(XT, K4); 
       XT:= "IF" LAST "THEN" XE "ELSE" X + H; HT:= 2 * H; 
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K1[J]:= ((((2 * K4[J] +  
       K2[J]) * .412022659166595 + K1[J]) * 2.23606797749979 -      
       K0[J]) * .375 - K3[J]) * HT + Y[J];      
       DER(XT, K1); 
       REJECT:= "FALSE"; FHM:= 0;     
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"      
       "BEGIN" DISCR:= ABS((1.6 * K3[J] - K2[J] - K4[J]) * 5 +      
           K0[J] + K1[J]);  
          TOL:= ABS(K0[J]) * E1 + E2; 
          REJECT:= DISCR > TOL "OR" REJECT;     
          FH:= DISCR / TOL; "IF" FH > FHM "THEN" FHM:= FH 
       "END";     
       MU:= 1 / (1 + FHM) + .45; "IF" REJECT "THEN"       
       "BEGIN" DATA[5]:= DATA[5] + 1; "IF" ABSH <= HMIN "THEN"      
          "BEGIN" DATA[6]:= DATA[6] + 1; HL:= H; REJECT:= "FALSE";  
             FIRST:= "TRUE"; "GOTO" NEXT
          "END";  
          H:= MU * H; "GOTO" TEST     
       "END";     
       "IF" FIRST "THEN"    
       "BEGIN" FIRST:= "FALSE"; HL:= H; H:= MU * H; "GOTO" ACC      
       "END";     
       FH:= MU * H / HL + MU - MU1; HL:= H; H:= FH * H;   
    ACC: MU1:= MU; HT:= HL / 12;      
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" Y[J]:= 
       ((K2[J] + K4[J]) * 5 + K0[J] + K1[J]) * HT + Y[J]; 
    NEXT: DATA[3]:= HL; DATA[4]:= DATA[4] + 1; X:= XT; OUT; 
       "IF" X ^= XE "THEN" "GOTO" TEST
    "END" RKE;    
         "EOP"    
1SECTION : 5.2.1.1.1.1.C      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : RK4A.
        
        
 AUTHOR:J.A.ZONNEVELD.      
        
        
 CONTRIBUTORS:M.BAKKER AND I.BRINK.   
        
        
 INSTITUTE: MATHEMATICAL CENTRE,      
        
        
 RECEIVED: 730715.
        
        
 BRIEF DESCRIPTION: 
        
     RK4A SOLVES AN INITIAL VALUE PROBLEM FOR A SINGLE FIRST ORDER  
     ORDINARY DIFFERENTIAL EQUATION   DY / DX = F(X,Y), WHERE F(X,Y)
     MAY BECOME LARGE, E.G. IN THE NEIGHBOURHOOD OF A SINGULARITY.  
     RK4A INTERCHANGES THE DEPENDENT AND INDEPENDENT VARIABLE.      
     THE EQUATION IS UPPOSED TO BE NON-STIFF.   
        
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEM, 
     SINGLE FIRST ORDER ORDINARY DIFFERENTIAL EQUATION    
     LARGE DERIVATIVE VALUES. 
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK4A(X,XA,B,Y,YA,FXY,E,D,FI,XDIR,POS);   
     "VALUE" FI,XDIR,POS;   
     "BOOLEAN" FI,XDIR,POS; 
     "REAL" X,XA,B,Y,YA,FXY;
     "ARRAY" E,D; 
     "CODE" 33016;
        
        
1SECTION : 5.2.1.1.1.1.C      (FEBRUARY 1979)                     PAGE 2      
        
        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
        
     X:    <VARIABLE>;      
           THE INDEPENDENT VARIABLE ; 
           UPON COMPLETION OF A CALL X IS EQUAL TO THE MOST RECENT  
           VALUE OF THE INDEPENDENT VARIABLE;   
     XA:   <ARITHMETIC EXPRESSION>;   
           ENTRY: THE INITIAL VALUE OF X;       
     B:    <ARITHMETIC EXPRESSION>;   
           B DEPENDS ON X AND Y.      
           THE EQUATION B=0 FULFILLED WITHIN THE TOLERANCES E[4] AND
           E[5]  SPECIFIES THE END OF THE INTEGRATION INTERVAL.     
           AT THE END OF EACH INTEGRATION STEP B IS EVALUATED AND IS
           TESTED FOR CHANGE OF SIGN; 
     Y:    <VARIABLE>;      
           THE DEPENDENT VARIABLE;    
     YA:   <ARITHMETIC EXPRESSION>;   
           ENTRY: THE VALUE OF Y AT X=XA;       
     FXY:  <ARITHMETIC EXPRESSION>;   
           FXY GIVES THE VALUE OF DY/DX;
           FXY USES X AND Y AS JENSEN PARAMETERS. 
     E:    <ARRAY IDENTIFIER>;
           "ARRAY" E[0:5];  
           ENTRY: 
           E[0], E[2] : RELATIVE TOLERANCES, FOR X AND Y RESPECTIVELY;
           E[1], E[3] : ABSOLUTE TOLERANCES, FOR X AND Y RESPECTIVELY;
           E[4], E[5] : TOLERANCES USED IN THE DETERMINATION OF     
           THE ZERO OF B;   
     D:    <ARRAY IDENTIFIER>;
           "ARRAY" D[0:4];  
           AFTER COMPLETION OF EACH STEP WE HAVE :
           IF D[0]>0 THEN X IS THE INTEGRATION VARIABLE;  
           IF D[0]<0 THEN Y IS THE INTEGRATION VARIABLE;  
           D[1] IS THE NUMBER OF STEPS SKIPPED; 
           D[2] IS THE STEPSIZE;      
           D[3] IS EQUAL TO THE LAST VALUE OF X;
           D[4] IS EQUAL TO THE LAST VALUE OF Y;
     FI:   <BOOLEAN EXPRESSION>;      
           IF FI="TRUE" THEN THE INTEGRATION IS STARTED WITH INITIAL
           VALUES X=XA,Y=YA.
           IF FI="FALSE" THEN THE INTEGRATION IS STARTED WITH X=D[3], 
           Y=D[4];
     XDIR, POS : <BOOLEAN EXPRESSION>;
           IF FI="TRUE" THEN THE INTEGRATION STARTS IN SUCH A WAY THAT
           IF POS="TRUE" AND XDIR="TRUE" THEN X INCREASES,
           IF POS="TRUE" AND XDIR="FALSE" THEN Y INCREASES, 
           IF POS="FALSE" AND XDIR="TRUE" THEN X DECREASES, 
           IF POS="FALSE" AND XDIR="FALSE" THEN Y DECREASES.
        
        
1SECTION : 5.2.1.1.1.1.C      (FEBRUARY 1979)                     PAGE 3      
        
        
        
 PROCEDURES USED : ZEROIN = CP34150.  
        
        
 METHOD AND PERFORMANCE:    
     RK4A IS BASED ON AN EXPLICIT, 5-TH ORDER RUNGE-KUTTA METHOD    
     AND INTERCHANGES THE DEPENDENT AND INDEPENDENT INTEGRATION VARIABLE.     
     IF   ABS(F(X,Y)) < 1, X IS USED AS INTEGRATION VARIABLE,OTHERWISE Y.     
     THE PROCEDURE IS PROVIDED WITH STEP SIZE AND ERROR CONTROL.    
     FOR DETAILS, E.G. HOW TO USE ARRAY E AND HOW TO SPECIFY THE ENDPOINT     
     SEE [1] ( RK4A IS A SLIGHTLY ADAPTED VERSION OF RK4).
        
        
 REFERENCES:      
     [1]J.A.ZONNEVELD,      
        AUTOMATIC NUMERICAL INTEGRATION.
        MATHEMATICAL CENTRE TRACT 8(1970).      
        
        
 EXAMPLE OF USE:  
        
     THE SOLUTION OF THE DIFFERENTIAL EQUATION  
     DY/DX=1-2*(X**2+Y), X>=0,
     Y=0               , X =0,
     IS REPRESENTED BY THE PARABOLA Y=X*(1-X);  
     WE WOULD LIKE TO FIND THE VALUE OF X FOR WHICH THE CURVE       
     OF THE SOLUTION INTERSECTS THE LINE Y+X=0. 
     THE SOLUTION CAN BE OBTAINED BY THE FOLLOWING PROGRAM: 
        
     "BEGIN" "COMMENT" INTEGRATION OF DY/DX=1-2*(Y+X**2), X>=0,     
         Y(0)=0,UNTIL THE CONDITION Y+X=0 IS SATISFIED;   
         "PROCEDURE" RK4A(X,XA,B,Y,YA,FXY,E,D,FI,XDIR,POS); 
         "CODE" 33016;      
         "REAL" X,Y;"ARRAY" D[0:4],E[0:5];      
         E[0]:=E[1]:=E[2]:=E[3]:=E[4]:=E[5]:="-6; 
         RK4A(X,0,X+Y,Y,0,1-2*(X*X+Y),E,D,"TRUE","TRUE","TRUE");    
         OUTPUT(61,"("10B"("X=")"+D.9D,10B"("EXACTLY:")"2B+D.9D/,   
         10B"("Y=")"+D.9D/,10B"("Y-X*(1-X)=")"+.10D")",X,2, 
         Y,Y-X*(1-X))       
     "END"
         "EOP"    
        
     THE PROGRAM PRINTS THE FOLLOWING RESULTS:  
        
     X=1.9999998554     EXACTLY: 2.0000000000   
        
     Y=-1.9999995347427     
        
     Y-X*(1-X)=0.0000000313 
        
1SECTION : 5.2.1.1.1.1.C      (AUGUST 1974)                       PAGE 4      
        
        
        
 SOURCE TEXT(S):  
0"CODE"" 33016 ;  
    "PROCEDURE" RK4A(X, XA, B, Y, YA, FXY, E, D, FI, XDIR,
    POS); "VALUE" FI, XDIR, POS; "BOOLEAN" FI, XDIR, POS; 
    "REAL" X, XA, B, Y, YA, FXY; "ARRAY" E, D;  
    "BEGIN" "INTEGER" I;    
       "BOOLEAN" IV, FIRST, FIR, REJ; 
       "REAL" K0, K1, K2, K3, K4, K5, FHM, ABSH, DISCR, S, XL,      
       COND0, S1, COND1, YL, HMIN, H, ZL, TOL, HL, MU, MU1; 
       "ARRAY" E1[1:2];     
        
       "BOOLEAN" "PROCEDURE" ZEROIN(X,Y,FX,EPS) ; "REAL" X,Y,FX,EPS ; 
       "CODE"            34150 ;      
        
       "PROCEDURE" RKSTEP(X, XL, H, Y, YL, ZL, FXY, D);   
       "VALUE" XL, YL, ZL, H; "REAL" X, XL, H, Y, YL, ZL, FXY;      
       "INTEGER" D; 
       "BEGIN" "IF" D = 2 "THEN" "GOTO" INTEGRATE; "IF" D = 3 "THEN"
          "BEGIN" X:= XL; Y:= YL; K0:= FXY * H "END"      
          "ELSE" "IF" D = 1 "THEN" K0:= ZL * H "ELSE" K0:= K0 * MU; 
          X:= XL + H / 4.5; Y:= YL + K0 / 4.5; K1:= FXY * H;
          X:= XL + H / 3; Y:= YL + (K0 + K1 * 3) / 12;    
          K2:= FXY * H; X:= XL + H * .5;
          Y:= YL + (K0 + K2 * 3) / 8; K3:= H * FXY;       
          X:= XL + H * .8;  
          Y:= YL + (K0 * 53 - K1 * 135 + K2 * 126 + K3 *  
          56) / 125; K4:= FXY * H; "IF" D <= 1 "THEN"     
          "BEGIN" X:= XL + H; 
             Y:= YL + (K0 * 133 - K1 * 378 + K2 * 276 + K3
             * 112 + K4 * 25) / 168; K5:= FXY * H;
             DISCR:= ABS(K0 * 21 - K2 * 162 + K3 * 224 - K4 
             * 125 + K5 * 42) / 14; "GOTO" END  
          "END";  
       INTEGRATE: X:= XL + H; 
          Y:= YL + ( - K0 * 63 + K1 * 189 - K2 * 36 - K3 *
          112 + K4 * 50) / 28; K5:= FXY * H;    
          Y:= YL + (K0 * 35 + K2 * 162 + K4 * 125 + K5 *  
          14) / 336;
       END:       
       "END" RKSTEP;                                           "COMMENT"      
1SECTION : 5.2.1.1.1.1.C      (AUGUST 1974)                       PAGE 5      
                                                                  ; 
        
        
       "REAL" "PROCEDURE" FZERO;      
       "BEGIN" "IF" IV "THEN" 
          "BEGIN" "IF" S = XL "THEN" FZERO:= COND0 "ELSE" "IF" S = S1 
             "THEN" FZERO:= COND1 "ELSE"
             "BEGIN" RKSTEP(X, XL, S - XL, Y, YL, ZL, FXY, 3);      
                FZERO:= B   
             "END"
          "END"   
          "ELSE"  
          "BEGIN" "IF" S = YL "THEN" FZERO:= COND0 "ELSE" "IF" S = S1 
             "THEN" FZERO:= COND1 "ELSE"
             "BEGIN" RKSTEP(Y, YL, S - YL, X, XL, ZL, 1 / 
                FXY, 3); FZERO:= B    
             "END"
          "END"   
       "END" FZERO; 
        
       "IF" FI "THEN"       
       "BEGIN" D[3]:= XA; D[4]:= YA; D[0]:= 1 "END";      
       D[1]:= 0; X:= XL:= D[3]; Y:= YL:= D[4]; IV:= D[0] > 0;       
       FIRST:= FIR:= "TRUE"; HMIN:= E[0] + E[1];
       H:= E[2] + E[3]; "IF" H < HMIN "THEN" HMIN:= H;    
    CHANGE: ZL:= FXY; "IF" ABS(ZL) <= 1 "THEN"  
       "BEGIN" "IF" "NOT"IV "THEN"    
          "BEGIN" D[2]:= H:= H / ZL; D[0]:= 1;  
             IV:= FIRST:= "TRUE"      
          "END";  
          "IF" FIR "THEN" "GOTO" A; I:= 1; "GOTO" AGAIN   
       "END"      
       "ELSE"     
       "BEGIN" "IF" IV "THEN" 
          "BEGIN" "IF" "NOT"FIR "THEN" D[2]:= H:= H * ZL; D[0]:= - 1; 
             IV:= "FALSE"; FIRST:= "TRUE"       
          "END";  
          "IF" FIR "THEN"   
          "BEGIN" H:= E[0] + E[1];    
          A: "IF" ("IF" FI "THEN" ("IF" IV "EQUIV" XDIR "THEN" H "ELSE"       
             H * ZL) < 0 "EQUIV" POS "ELSE" H * D[2] < 0) "THEN" H:= - H      
          "END";  
          I:= 1   
       "END";                                                  "COMMENT"      
1SECTION : 5.2.1.1.1.1.C      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
    AGAIN: ABSH:= ABS(H); "IF" ABSH < HMIN "THEN" 
       "BEGIN" H:= SIGN(H) * HMIN; ABSH:= HMIN "END";     
       "IF" IV "THEN"       
       "BEGIN" RKSTEP(X, XL, H, Y, YL, ZL, FXY, I);       
          TOL:= E[2] * ABS(K0) + E[3] * ABSH    
       "END"      
       "ELSE"     
       "BEGIN" RKSTEP(Y, YL, H, X, XL, 1 / ZL, 1 / FXY, I); 
          TOL:= E[0] * ABS(K0) + E[1] * ABSH    
       "END";     
       REJ:= DISCR > TOL; MU:= TOL / (TOL + DISCR) + .45; 
       "IF" REJ "THEN"      
       "BEGIN" "IF" ABSH <= HMIN "THEN" 
          "BEGIN" "IF" IV "THEN"      
             "BEGIN" X:= XL + H; Y:= YL + K0 "END"
             "ELSE" 
             "BEGIN" X:= XL + K0; Y:= YL + H "END";       
             D[1]:= D[1] + 1; FIRST:= "TRUE"; "GOTO" NEXT 
          "END";  
          H:= H * MU; I:= 0; "GOTO" AGAIN       
       "END" REJ; 
       "IF" FIRST "THEN"    
       "BEGIN" FIRST:= FIR; HL:= H; H:= MU * H; "GOTO" ACCEPT       
       "END";     
       FHM:= MU * H / HL + MU - MU1; HL:= H; H:= FHM * H; 
    ACCEPT: "IF" IV "THEN" RKSTEP(X, XL, HL, Y, YL, ZL, FXY,
       2) "ELSE" RKSTEP(Y, YL, HL, X, XL, ZL, 1 / FXY, 2);
       MU1:= MU;  
    NEXT: "IF" FIR "THEN"   
       "BEGIN" FIR:= "FALSE"; COND0:= B;
          "IF" "NOT"(FI "OR" REJ) "THEN" H:= D[2] 
       "END"      
       "ELSE"     
       "BEGIN" D[2]:= H; COND1:= B;   
          "IF" COND0 * COND1 <= 0 "THEN" "GOTO" ZERO;     
          COND0:= COND1     
       "END";     
       D[3]:= XL:= X; D[4]:= YL:= Y; "GOTO" CHANGE;       
    ZERO: E1[1]:= E[4]; E1[2]:= E[5]; 
       S1:= "IF" IV "THEN" X "ELSE" Y;
       S:= "IF" IV "THEN" XL "ELSE" YL ;
       ZEROIN(S,S1,FZERO,ABS(E1[1]*S)+ABS(E1[2])) ;       
       S1:= "IF" IV "THEN" X "ELSE" Y ; 
       "IF" IV "THEN" RKSTEP(X, XL, S - XL, Y, YL, ZL, FXY, 3)      
       "ELSE" RKSTEP(Y, YL, S - YL, X, XL, ZL, 1 / FXY,   
       3); D[3]:= X; D[4]:= Y 
    "END" RK4A;   
         "EOP"    
1SECTION : 5.2.1.1.1.1.D      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : RK4NA. 
        
        
 AUTHOR:J.A.ZONNEVELD.      
        
        
 CONTRIBUTORS:M.BAKKER AND I.BRINK.   
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730715.
        
        
 BRIEF DESCRIPTION: 
        
     RK4NA SOLVES AN INITIAL VALUE PROBLEM FOR A SYSTEM OF FIRST    
     ORDER ORDINARY DIFFERENTIAL EQUATIONS  DY / DX = F(X,Y), OF WHICH
     THE DERIVATIVE COMPONENTS ARE SUPPOSED TO BECOME LARGE, E.G. IN
     THE NEIGHBOURHOOD OF SINGULARITIES. RK4NA INTERCHANGES THE     
     INDEPENDENT VARIABLE AND ONE DEPENDENT VARIABLE. THE SYSTEM    
     IS SUPPOSED TO BE NON-STIFF.     
        
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEM, 
     SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS, 
     LARGE DERIVATIVE COMPONENTS.     
        
        
1SECTION : 5.2.1.1.1.1.D      (FEBRUARY 1979)                     PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK4NA(X,XA,B,FXJ,J,E,D,FI,N,L,POS);      
     "VALUE" FI,N,L,POS;    
     "INTEGER" J,N,L;       
     "BOOLEAN" FI,POS;      
     "REAL" B,FXJ;
     "ARRAY" X,XA,E,D;      
     "CODE" 33017;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
        
     X:    <ARRAY IDENTIFIER>;
           "ARRAY" X[0:N];  
           X[0] IS THE INDEPENDENT VARIABLE,    
           X[1],...,X[N] ARE THE DEPENDENT VARIABLES;     
           EXIT:THE SOLUTION AT B=0;  
     XA:   <ARRAY IDENTIFIER>;
           "ARRAY" XA[0:N]; 
           ENTRY: THE INITIAL VALUES OF X[0],...,X[N];    
     B:    <ARITHMETIC EXPRESSION>;   
           B DEPENDS ON X[0],...,X[N];
           IF THE EQUATION B=0 IS SATISFIED WITHIN A      
           CERTAIN TOLERANCE,THE INTEGRATION IS TERMINATED; 
           B IS EVALUATED AND TESTED FOR CHANGE OF SIGN AT THE      
           END OF EACH STEP;
           FOR THE TOLERANCE SEE PARAMETER E.   
     FXJ:  <ARITHMETIC EXPRESSION>;   
           FXJ DEPENDS ON X[0],...,X[N] AND J, DEFINING THE RIGHT   
           HAND SIDE OF THE DIFFERENTIAL EQUATION;
           AT EACH CALL IT DELIVERS :DX[J]/DX[0]; 
     J:    <VARIABLE>;      
           J IS USED AS A JENSEN PARAMETER FOR FXJ;       
     E:    <ARRAY IDENTIFIER>;
           "ARRAY" E[0:2*N+3];
           ENTRY: E[2*J] AND E[2*J+1] , 0<=J<=N , 
           ARE THE RELATIVE AND THE ABSOLUTE TOLERANCE ,  
           RESPECTIVELY, ASSOCIATED WITH X[J];  
           E[2*N+2] AND E[2*N+3] ARE THE RELATIVE AND ABSOLUTE      
           TOLERANCE  USED IN THE DETERMINATION OF THE ZERO OF B;   
     D:    <ARRAY IDENTIFIER>;
           "ARRAY" D[0:N+3];
           AFTER COMPLETION OF EACH STEP WE HAVE :
           ENTIER(D[0]+.5) IS THE NUMBER OF STEPS SKIPPED;
           D[2] IS THE STEP LENGTH;   
           D[J+3] IS THE LAST VALUE OF X[J], J=0,...,J=N; 
1SECTION : 5.2.1.1.1.1.D      (DECEMBER 1975)                     PAGE 3      
        
        
        
     FI:   <BOOLEAN EXPRESSION>       
           IF FI="TRUE" THEN THE INTEGRATION IS STARTED WITH INITIAL
           CONDITIONS X[J]=XA[J];IF FI="FALSE" THEN THE INTEGRATION IS
           CONTINUED WITH X[J]=D[J+3];
     N:    <ARITHMETIC EXPRESSION>;   
           THE NUMBER OF EQUATIONS;   
     L:    <ARITHMETIC EXPRESSION>;   
           AN INTEGER TO BE SUPPLIED BY THE USER, 0<=L<=N (SEE POS);
     POS:  <BOOLEAN EXPRESSION>       
           IF FI="TRUE" THEN THE INTEGRATION STARTS IN SUCH A WAY   
           THAT X[L] INCREASES IF POS="TRUE" AND X[L] DECREASES IF  
           POS="FALSE";     
           IF FI="FALSE" THEN POS IS OF NO SIGNIFICANCE.  
        
        
 PROCEDURES USED : ZEROIN = CP34150.  
        
        
 REQUIRED CENTRAL MEMORY : CIRCA 9*N MEMORY PLACES.       
        
        
 METHOD AND PERFORMANCE :   
     RK4NA IS BASED ON AN EXPLICIT, 5-TH ORDER RUNGE-KUTTA METHOD   
     AND INTERCHANGES VARIABLES. THE PROCEDURE IS PROVIDED WITH STEPSIZE      
     AND ERROR CONTROL. FOR DETAILS, E.G. HOW TO USE ARRAY E, HOW TO
     SPECIFY THE ENDPOINT, HOW TO USE L AND POS, SEE [1] ( RK4NA IS A 
     SLIGHTLY ADAPTED VERSION OF RK4N ).
        
        
 REFERENCES:      
        
     [1].J.A.ZONNEVELD.     
         AUTOMATIC NUMERICAL INTEGRATION.       
         MATHEMATICAL CENTRE TRACT 8 (1970).    
        
        
1SECTION : 5.2.1.1.1.1.D      (FEBRUARY 1979)                     PAGE 4      
        
        
        
 EXAMPLE OF USE:  
        
     THE PERIOD OF THE SOLUTION OF THE VAN DER POL EQUATION 
     DX[1]/DT = X[2]
     DX[2]/DT = 10*(1-X[1]**2)*X[2]-X[1], T>=0, 
     CAN BE OBTAINED BY THE FOLLOWING PROGRAM:  
        
     "BEGIN" "COMMENT" VAN DER POL;   
         "PROCEDURE" RK4NA(X,XA,B,FXY,J,E,D,FI,N,L,POS);  
         "CODE" 33017;      
         "REAL" X0; 
        
         "PROCEDURE" PRINT(X);"ARRAY" X;
         OUTPUT(61,"("/,4(+ZD.10D3B)")",X[0],X[1],X[2],X0); 
        
         "INTEGER" J,K;"BOOLEAN" FIRST; 
         "ARRAY" E[0:7],XA,X[0:2],D[0:5];       
         "FOR" K:=0,1,2,3,4,5 "DO" E[K]:=.1"-6; E[6]:=E[7]:="-8 ;   
         OUTPUT(61,"(""("VAN DER POL")",/BB"("EPS=")"D.10D,/BB      
         "("THE VALUES OF X[0],X[1],X[2],P,RESPECTIVELY")"")",E[0]);
         X0:=XA[0]:=XA[2]:=0;XA[1]:=2;J:=0;PRINT(XA);     
         FIRST:="TRUE";     
         "FOR" J:=1,2,3,4 "DO"
         "BEGIN" RK4NA(X,XA,X[2],"IF" K=1 "THEN" X[2] "ELSE"
             10*(1-X[1]**2)*X[2]-X[1],K,E,D,FIRST,2,0,"TRUE");      
             X0:=X[0]-X0;PRINT(X);FIRST:="FALSE"; X0:=X[0]
         "END"    
     "END"
         "EOP"    
        
        
     THE PROGRAM PRINTS THE FOLLOWING RESULTS:  
        
     VAN DER POL  
        
     EPS=0.0000001000       
        
     THE VALUES OF X[0],X[1],X[2],P,RESPECTIVELY: 
        
     +0.00000000  +2.00000000  +0.00000000  +0.00000000   
        
     +9.32386570  -2.01428560  +0.00000000  +9.32386570   
        
     +18.86305411 +2.01428557  +0.00000000  +9.53918840   
        
     +28.40224194  -2.01428858  -0.00000000  +9.53918783  
        
     +37.94143003  +2.01428558  +0.00000000  +9.53918809  
        
1SECTION : 5.2.1.1.1.1.D      (AUGUST 1974)                       PAGE 5      
        
        
        
 SOURCE TEXT(S):  
0"CODE"" 33017 ;  
    "PROCEDURE" RK4NA(X, XA, B, FXJ, J, E, D, FI, N, L, POS);       
    "VALUE" FI, N, L, POS; "INTEGER" J, N, L; "BOOLEAN" FI, POS;    
    "REAL" B, FXJ; "ARRAY" X, XA, E, D; 
    "BEGIN" "INTEGER" I, IV, IV0;     
       "BOOLEAN" FIR, FIRST, REJ;     
       "REAL" H, COND0, COND1, FHM, ABSH, TOL, FH, MAX, X0, 
       X1, S, HMIN, HL, MU, MU1;      
       "ARRAY" XL, DISCR, Y[0:N], K[0:5,0:N], E1[1:2];    
        
       "BOOLEAN" "PROCEDURE" ZEROIN(X,Y,FX,EPS) ; "REAL" X,Y,FX,EPS ; 
       "CODE"            34150 ;      
        
       "PROCEDURE" RKSTEP(H, D); "VALUE" H, D; "INTEGER" D; "REAL" H; 
       "BEGIN" "INTEGER" I; 
        
        
          "PROCEDURE" F(T); "VALUE" T; "INTEGER" T;       
          "BEGIN" "INTEGER" I;
             "REAL" P;      
             "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" Y[J]:= FXJ;
             P:= H / Y[IV]; 
             "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" "IF" I ^= IV "THEN" K[T,I]:= Y[I] * P "END"    
          "END" F;
        
          "IF" D = 2 "THEN" "GOTO" INTEGRATE; "IF" D = 3 "THEN"     
          "BEGIN" "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I]; 
             F(0) 
          "END"   
          "ELSE" "IF" D = 1 "THEN"    
          "BEGIN" "REAL" P; 
             P:= H / Y[IV]; 
             "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" "IF" I ^= IV "THEN" K[0,I]:= P * Y[I] "END"    
          "END"   
          "ELSE"  
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"   
          "BEGIN" "IF" I ^= IV "THEN" K[0,I]:= K[0,I] * MU "END";   
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + ("IF" I
          = IV "THEN" H "ELSE" K[0,I]) / 4.5; F(1);       
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + ("IF" I
          = IV "THEN" H * 4 "ELSE" (K[0,I] + K[1,I] * 3)) / 12;     
          F(2);   
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + ("IF" I
          = IV "THEN" H * .5 "ELSE" (K[0,I] + K[2,I] * 3) / 8);     
          F(3);                                                "COMMENT"      
1SECTION : 5.2.1.1.1.1.D      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + ("IF" I
          = IV "THEN" H * .8 "ELSE" (K[0,I] * 53 - K[1,I] * 135     
          + K[2,I] * 126 + K[3,I] * 56) / 125); F(4);     
          "IF" D <= 1 "THEN"
          "BEGIN" "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] +
             ("IF" I = IV "THEN" H "ELSE" (K[0,I] * 133 - 
             K[1,I] * 378 + K[2,I] * 276 + K[3,I] * 112 + 
             K[4,I] * 25) / 168); F(5); 
             "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" "IF" I ^= IV "THEN" DISCR[I]:= ABS(K[0,I] * 21 
                - K[2,I] * 162 + K[3,I] * 224 - K[4,I] *  
                125 + K[5,I] * 42) / 14 
             "END"; 
             "GOTO" END     
          "END";  
       INTEGRATE: "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I]  
          + ("IF" I = IV "THEN" H "ELSE" ( - K[0,I] * 63 + K[1,I]   
          * 189 - K[2,I] * 36 - K[3,I] * 112 + K[4,I] * 50) 
          / 28); F(5);      
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"   
          "BEGIN" "IF" I ^= IV "THEN" X[I]:= XL[I] + (K[0,I] * 35   
             + K[2,I] * 162 + K[4,I] * 125 + K[5,I] * 14) / 336     
          "END" ; 
          END ..  
          "END" RKSTEP ;    
        
        
       "REAL" "PROCEDURE" FZERO;      
       "BEGIN" "IF" S = X0 "THEN" FZERO:= COND0 "ELSE" "IF" S = X1  
          "THEN" FZERO:= COND1 "ELSE" 
          "BEGIN" RKSTEP(S - XL[IV], 3); FZERO:= B "END"  
       "END" FZERO; 
        
       "IF" FI "THEN"       
       "BEGIN" "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" D[I + 3]:= XA[I];
          D[0]:= D[2]:= 0   
       "END";     
       D[1]:= 0;  
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I]:= D[I + 3]; 
       IV:= D[0]; H:= D[2]; FIRST:= FIR:= "TRUE"; Y[0]:= 1; 
       "GOTO" CHANGE;       
    AGAIN: ABSH:= ABS(H); "IF" ABSH < HMIN "THEN" 
       "BEGIN" H:= "IF" H > 0 "THEN" HMIN "ELSE" - HMIN;  
          ABSH:= ABS(H)     
       "END";     
       RKSTEP(H, I); REJ:= "FALSE"; FHM:= 0;    
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"      
       "BEGIN" "IF" I ^= IV "THEN"    
          "BEGIN" TOL:= E[2 * I] * ABS(K[0,I]) + E[2 * I + 1]       
             * ABSH; REJ:= TOL < DISCR[I] "OR" REJ;       
             FH:= DISCR[I] / TOL; "IF" FH > FHM "THEN" FHM:= FH     
          "END"   
       "END";                                                  "COMMENT"      
1SECTION : 5.2.1.1.1.1.D      (AUGUST 1974)                       PAGE 7      
                                                                  ; 
        
        
       MU:= 1 / (1 + FHM) + .45; "IF" REJ "THEN"
       "BEGIN" "IF" ABSH <= HMIN "THEN" 
          "BEGIN" "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"     
             "BEGIN" "IF" I ^= IV "THEN" X[I]:= XL[I] + K[0,I]      
                "ELSE" X[I]:= XL[I] + H 
             "END"; 
             D[1]:= D[1] + 1; FIRST:= "TRUE"; "GOTO" NEXT 
          "END";  
          H:= H * MU; I:= 0; "GOTO" AGAIN       
       "END";     
       "IF" FIRST "THEN"    
       "BEGIN" FIRST:= FIR; HL:= H; H:= MU * H; "GOTO" ACCEPT       
       "END";     
       FH:= MU * H / HL + MU - MU1; HL:= H; H:= FH * H;   
    ACCEPT: RKSTEP(HL, 2); MU1:= MU;  
    NEXT: "IF" FIR "THEN"   
       "BEGIN" FIR:= "FALSE"; COND0:= B;
          "IF" "NOT"(FI "OR" REJ) "THEN" H:= D[2] 
       "END"      
       "ELSE"     
       "BEGIN" D[2]:= H; COND1:= B;   
          "IF" COND0 * COND1 <= 0 "THEN" "GOTO" ZERO;     
          COND0:= COND1     
       "END";     
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" D[I + 3]:= XL[I]:= X[I]; 
    CHANGE: IV0:= IV;       
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" Y[J]:= FXJ;    
       MAX:= ABS(Y[IV]);    
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"      
       "BEGIN" "IF" ABS(Y[I]) > MAX "THEN"      
          "BEGIN" MAX:= ABS(Y[I]); IV:= I "END" 
       "END";     
       "IF" IV0 ^= IV "THEN"
       "BEGIN" FIRST:= "TRUE"; D[0]:= IV;       
          D[2]:= H:= Y[IV] / Y[IV0] * H 
       "END";     
       X0:= XL[IV]; "IF" FIR "THEN"   
       "BEGIN" HMIN:= E[0] + E[1];    
          "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"   
          "BEGIN" H:= E[2 * I] + E[2 * I + 1];  
             "IF" H < HMIN "THEN" HMIN:= H      
          "END";  
          H:= E[2 * IV] + E[2 * IV + 1];
             "IF" (FI "AND" (Y[L]/Y[IV]*H<0 "EQUIV" POS)) "OR"      
             ("NOT" FI "AND" D[2]*H<0) "THEN" H:= -H      
       "END";     
       I:= 1; "GOTO" AGAIN; 
    ZERO: E1[1]:= E[2 * N + 2]; E1[2]:= E[2 * N + 3];     
       X1:=X[IV] ; S:=X0 ;  
       ZEROIN(S,X1,FZERO,ABS(E1[1]*S) + ABS(E1[2])) ; X0:=S ; X1:=X[IV];      
       RKSTEP(X0 - XL[IV], 3);
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" D[I + 3]:= X[I]
    "END" RK4NA;  
         "EOP"    
1SECTION : 5.2.1.1.1.1.E      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : RK5NA. 
        
        
 AUTHOR: J.A.ZONNEVELD.     
        
        
 CONTRIBUTORS: M.BAKKER AND I.BRINK.  
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730715.
        
        
 BRIEF DESCRIPTION: 
        
     RK5NA SOLVES AN INITIAL VALUE PROBLEM FOR A SYSTEM OF FIRST ORDER
     ORDINARY DIFFERENTIAL EQUATIONS  DY / DX = F(X,Y), OF WHICH THE
     DERIVATIVE COMPONENTS ARE SUPPOSED TO BECOME LARGE, E.G. IN THE
     NEIGHBOURHOOD OF SINGULARITIES. RK5NA INTRODUCES THE ARC LENGTH
     AS INTEGRATION VARIABLE. THE SYSTEM IS SUPPOSED TO BE NON-STIFF. 
        
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEM. 
     SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS. 
     LARGE DERIVATIVE COMPONENTS.     
        
        
1SECTION : 5.2.1.1.1.1.E      (FEBRUARY 1979)                     PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK5NA(X, XA, B, FXJ, J, E, D, FI, N, L, POS);      
     "VALUE" FI, N, L, POS; "INTEGER" J, N, L; "BOOLEAN" FI, POS;   
     "REAL" B, FXJ; "ARRAY" X, XA, E, D;
     "CODE" 33018;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
        
     X:  <ARRAY IDENTIFIER>;
         "ARRAY" X[0 : N];  
         THE DEPENDENT VARIABLES;     
         X[0], ..., X[N] CAN BE USED AS JENSEN PARAMETERS;
     XA:  <ARRAY IDENTIFIER>; 
         "ARRAY" XA[0 : N]; 
         ENTRY: THE INITIAL VALUES OF X[J], J = 0, ..., N;
     B:  <ARITHMETIC EXPRESSION>;     
         B DEPENDS ON X[0], ..., X[N];
         IF,WITHIN SOME TOLERANCE,B = 0 THEN THE INTEGRATION IS TER-
         MINATED; SEE ALSO THE EXPLANATION OF THE PARAMETER E;      
     FXJ:  <ARITHMETIC EXPRESSION>;   
         THE RIGHT HAND SIDE OF THE DIFFERENTIAL EQUATION;
         FXJ DEPENDS ON X[0], ..., X[N], J, GIVING THE VALUE OF     
         DX[J] / DX[0];     
     J:  <VARIABLE>;
         J IS USED AS A JENSEN PARAMETER TO DENOTE,IN THE ACTUAL    
         PARAMETER CORRESPONDING TO FXJ,THE NUMBER OF THE FUNCTION  
         REQUIRED;
     E:  <ARRAY IDENTIFIER>;
         "ARRAY" E[0 : 2 * N + 3];    
         ENTRY:   
         E[2 * J] AND E[2 * J + 1] ARE THE RELATIVE AND THE ABSOLUTE
         TOLERANCE,RESPECTIVELY,ASSOCIATED WITH X[J],J = 0,..., N, WHILE      
         E[2 * N + 2] AND E[2 * N + 3] ARE THE ONES ASSOCIATED WITH B;
     D:  <ARRAY IDENTIFIER>;
         "ARRAY" D[1 : N + 3];
         AFTER COMPLETION OF EACH STEP WE HAVE: 
        
             ABS(D[1])                   THE ARC LENGTH,  
                 D[2]                    THE STEP LENGTH, 
                 D[I + 3]                THE LATEST VALUE OF X[I],  
                                         I = 0, ..., N;   
        
     FI:  <BOOLEAN EXPRESSION>;       
         IF FI = "TRUE" THEN THE INTEGRATION IS STARTED WITH INITIAL
         CONDITIONS X[I] = XA[I], I = 0, ..., N;
         IF FI = "FALSE" THEN THE INTEGRATION IS CONTINUED WITH     
         X[I] = D[I + 3];   
     N:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF EQUATIONS;     
1SECTION : 5.2.1.1.1.1.E      (FEBRUARY 1979)                     PAGE 3      
        
        
        
     L:  <ARITHMETIC EXPRESSION>;     
         1 <= L <= N; SEE THE EXPLANATION OF POS; 
     POS:  <BOOLEAN EXPRESSION>;      
         IF FI = "TRUE" THEN THE INTEGRATION STARTS IN SUCH A WAY THAT
         X[L] INCREASES IF POS = "TRUE" AND DECREASES IF POS = "FALSE".       
         IF FI = "FALSE" THEN POS IS IGNORED.   
        
        
 PROCEDURES USED : ZEROIN = CP34150.  
        
        
 REQUIRED CENTRAL MEMORY : CIRCA 9 * N  MEMORY PLACES.    
        
        
 METHOD AND PERFORMANCE :   
     RK5NA IS BASED ON A 5-TH ORDER RUNGE-KUTTA  METHOD AND INTEGRATES
     THE SYSTEM  DX[J] / DX[0] = F(J,X[0],..,X[N]) / F(0,X[0],..,X[N]).       
     THE ARC LENGTH IS INTRODUCED AS INTEGRATION VARIABLE.
     THE INTEGRATION PROCESS IS TERMINATED IF SOME CONDITION ON     
     X[0], .. ,X[N], TO BE SUPPLIED BY THE USER, IS SATISFIED.      
     RK5NA USES STEPLENGTH AND ERROR CONTROL. DETAILS ABOUT THE     
     PROCEDURE AND THE UNDERLYING THEORY ARE GIVEN IN [1] ( RK5NA IS
     A SLIGHTLY  ADAPTED VERSION OF RK5N).      
        
        
 REFERENCES:      
        
     [1]. J.A.ZONNEVELD,    
          AUTOMATIC NUMERICAL INTEGRATION,      
          MATH.CENTRE TRACT 8 (1970) .
        
        
1SECTION : 5.2.1.1.1.1.E      (AUGUST 1974)                       PAGE 4      
        
        
        
 EXAMPLE OF USE:  
        
     THE VAN DER POL EQUATION IN THE PHASE PLANE
        
         DX[1] / DX[0] = (10*(1-X[0]**2)*X[1]-X[0])/X[1]  
        
     CAN BE INTEGRATED BY THE PROCEDURE RK5NA; THE STARTING VALUES ARE
     X[0] = 2, X[1] = 0. THE INTEGRATION PROCEEDS UNTIL THE NEXT ZERO OF      
     X[1],THEN IT CONTINUES UNTIL THE NEXT ZERO AND SO ON UNTIL THE 
     FOURTH ZERO IS ENCOUNTERED.IN THE EXAMPLE THE OUTPUT IS GIVEN FOR
     THE TOLERANCES E[0]=E[1]=E[2]=E[3]= "-6, E[4]=E[5]= "-10.      
     THE PROGRAM READS AS FOLLOWS:    
        
     "BEGIN" "COMMENT" VAN DER POL IN THE PHASE PLANE;    
         "PROCEDURE" RK5NA(X,XA,B,FXJ,J,E,D,FI,N,L,POS); "CODE" 33018;
         "INTEGER" J,K; "BOOLEAN" FIRST;
         "ARRAY" E[0:5],XA,X[0:1],D[1:4];       
         "PROCEDURE" PRINT(X); "ARRAY" X;       
         "BEGIN" OUTPUT(61,"("/B"("X[0]=")"+D.10D,
             10B"("X[1]=")"+D.10D,10B"("S=")"3D.10D")",X[0],
             X[1],ABS(D[1]))
         "END";   
         "FOR" K:=0,1,2,3 "DO" E[K]:="-6 ; E[4]:=E[5]:="-10 ;  D[1]:=0;       
         XA[0]:=2; XA[1]:=0; J:=0; PRINT(XA); AA: 
         FIRST:=J=0;
         RK5NA(X,XA,X[1],"IF" K=0 "THEN" X[1] "ELSE"      
         10*(1-X[0]**2)*X[1]-X[0],K,E,D,FIRST,1,1,"FALSE");;J:=J+1; 
         PRINT(X); "IF" J<4 "THEN" "GO TO" AA   
     "END"
         "EOP"    
        
     RESULTS:     
        
     X[0]=+2.0000000000    X[1]=+0.0000000000    S=000.0000000000   
        
     X[0]=-2.0142853657    X[1]=-0.0000000012    S=029.3873834087   
        
     X[0]=+2.0142853659    X[1]=+0.0000000001    S=058.7884331939   
        
     X[0]=-2.0142853659    X[1]=-0.0000000000    S=088.1894829781   
        
     X[0]=+2.0142853659    X[1]=+0.0000000000    S=117.5905327623   
        
1SECTION : 5.2.1.1.1.1.E      (AUGUST 1974)                       PAGE 5      
        
        
        
 SOURCE TEXT(S):  
 "CODE"" 33018 ;  
    "PROCEDURE" RK5NA(X, XA, B, FXJ, J, E, D, FI, N, L, POS);       
    "VALUE" FI, N, L, POS; "INTEGER" J, N, L; "BOOLEAN" FI, POS;    
    "REAL" B, FXJ; "ARRAY" X, XA, E, D; 
    "BEGIN" "INTEGER" I;    
       "BOOLEAN" FIRST, FIR, REJ;     
       "REAL" FHM, S, S0, COND0, S1, COND1, H, ABSH, TOL, FH,       
       HL, MU, MU1; 
       "ARRAY" Y, XL, DISCR[0:N], K[0:5,0:N], E1[1:2];    
       "REAL" "PROCEDURE" SUM(J,A,B,XJ) ; "INTEGER" J,A,B ; "REAL" XJ ;       
       "BEGIN" "REAL" S ; S:= 0 ;     
          "FOR" J:=A "STEP" 1 "UNTIL" B "DO" S:=S+XJ ; SUM:= S      
       "END" SUM ;
       "BOOLEAN" "PROCEDURE" ZEROIN(X,Y,FX,EPS) ; "REAL" X,Y,FX,EPS ; 
       "CODE"            34150 ;      
       "PROCEDURE" RKSTEP(H, D); "VALUE" H, D; "INTEGER" D; "REAL" H; 
       "BEGIN" "INTEGER" I; 
          "PROCEDURE" F(T); "VALUE" T; "INTEGER" T;       
          "BEGIN" "INTEGER" I;
             "REAL" P;      
             "FOR" J:= 0 "STEP" 1 "UNTIL" N "DO" Y[J]:= FXJ;
             P:= H / SQRT(SUM(I, 0, N, Y[I] ** 2));       
             "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" K[T,I]:= Y[I] * P  
          "END" F;
          "IF" D = 2 "THEN" "GOTO" INTEGRATE; "IF" D = 1 "THEN"     
          "BEGIN" "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" K[0,I]:= K[0,I] 
             * MU; "GOTO" A 
          "END";  
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I]; F(0);   
       A: "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] +
          K[0,I] / 4.5; F(1); 
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + (K[0,I]
          + K[1,I] * 3) / 12; F(2);   
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + (K[0,I]
          + K[2,I] * 3) / 8; F(3);    
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + (K[0,I]
          * 53 - K[1,I] * 135 + K[2,I] * 126 + K[3,I] * 56) 
          / 125; F(4); "IF" D <= 1 "THEN"       
          "BEGIN" "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] +
             (K[0,I] * 133 - K[1,I] * 378 + K[2,I] * 276 +
             K[3,I] * 112 + K[4,I] * 25) / 168; F(5);     
             "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" DISCR[I]:= 
             ABS(K[0,I] * 21 - K[2,I] * 162 + K[3,I] * 224
             - K[4,I] * 125 + K[5,I] * 42) / 14; "GOTO" END 
          "END";  
       INTEGRATE: "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I]  
          + ( - K[0,I] * 63 + K[1,I] * 189 - K[2,I] * 36 -
          K[3,I] * 112 + K[4,I] * 50) / 28; F(5); 
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I] + (K[0,I]
          * 35 + K[2,I] * 162 + K[4,I] * 125 + K[5,I] * 14) 
          / 336;  
       END:       
       "END" RKSTEP;                                           "COMMENT"      
1SECTION : 5.2.1.1.1.1.E      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
       "REAL" "PROCEDURE" FZERO;      
       "BEGIN" "IF" S = S0 "THEN" FZERO:= COND0 "ELSE" "IF" S = S1  
          "THEN" FZERO:= COND1 "ELSE" 
          "BEGIN" RKSTEP(S - S0, 3); FZERO:= B "END"      
       "END" FZERO; 
        
       "IF" FI "THEN"       
       "BEGIN" "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" D[I + 3]:= XA[I];
          D[1]:= D[2]:= 0   
       "END";     
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= XL[I]:= D[I + 3]; 
       S:= D[1]; FIRST:= FIR:= "TRUE"; H:= E[0] + E[1];   
       "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"      
       "BEGIN" ABSH:= E[2 * I] + E[2 * I + 1];  
          "IF" H > ABSH "THEN" H:= ABSH 
       "END";     
       "IF" FI "THEN"       
       "BEGIN" J:= L; "IF" FXJ * H < 0 "EQUIV" POS "THEN" H:= - H "END"       
       "ELSE" "IF" D[2] * H < 0 "THEN" H:= - H; I:= 0;    
    AGAIN: RKSTEP(H, I); REJ:= "FALSE"; FHM:= 0;
       ABSH:= ABS(H);       
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"      
       "BEGIN" TOL:= E[2 * I] * ABS(K[0,I]) + E[2 * I + 1] *
          ABSH; REJ:= TOL < DISCR[I] "OR" REJ;  
          FH:= DISCR[I] / TOL; "IF" FH > FHM "THEN" FHM:= FH
       "END";     
       MU:= 1 / (1 + FHM) + .45; "IF" REJ "THEN"
       "BEGIN" H:= H * MU; I:= 1; "GOTO" AGAIN "END";     
       "IF" FIRST "THEN"    
       "BEGIN" FIRST:= FIR; HL:= H; H:= MU * H "END"      
       "ELSE"     
       "BEGIN" FH:= MU * H / HL + MU - MU1; HL:= H; H:= FH * H      
       "END";     
    ACCEPT: RKSTEP(HL, 2); MU1:= MU; S:= S + HL;
       "IF" FIR "THEN"      
       "BEGIN" COND0:= B; FIR:= "FALSE"; "IF" "NOT"FI "THEN" H:= D[2] 
       "END"      
       "ELSE"     
       "BEGIN" D[2]:= H; COND1:= B;   
          "IF" COND0 * COND1 <= 0 "THEN" "GOTO" ZERO;     
          COND0:= COND1     
       "END";     
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" D[I + 3]:= XL[I]:= X[I]; 
       D[1]:= S0:= S; I:= 0; "GOTO" AGAIN;      
    ZERO: E1[1]:= E[2 * N + 2]; E1[2]:= E[2 * N + 3];     
       S1:=S ; S:=S0 ;      
       ZEROIN(S,S1,FZERO,ABS(E1[1]*S)+ABS(E1[2])) ;       
       RKSTEP(S - S0, 3);   
       "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" D[I + 3]:= X[I]; D[1]:= S
    "END" RK5NA;  
         "EOP"    
1SECTION : 5.2.1.1.1.1.F      (FEBRUARY 1979)                     PAGE 1
 
 
 
 PROCEDURE : MULTISTEP.
 
 
 AUTHOR: P.W.HEMKER.
 
 
 INSTITUTE: MATHEMATICAL CENTRE.
 
 
 RECEIVED: 730515.
 
 
 BRIEF DESCRIPTION:
 
     MULTISTEP  SOLVES  AN INITIAL VALUE  PROBLEM, FOR A SYSTEM OF
     FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS  DY = F(Y).
     IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE INTEGRATION OF
     STIFF DIFFERENTIAL EQUATIONS. IT CAN ALSO BE USED FOR
     NON-STIFF PROBLEMS.
 
 
 KEYWORDS:
 
     INITIAL VALUE PROBLEM,
     SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS.
     STIFF EQUATIONS,
 
 
1SECTION : 5.2.1.1.1.1.F      (FEBRUARY 1979)                     PAGE 2
 
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "BOOLEAN" "PROCEDURE" MULTISTEP(X,XEND,Y,HMIN,HMAX,YMAX,EPS,
                       FIRST,SAVE,DERIV,AVAILABLE,JACOBIAN,STIFF,N,OUT);
     "VALUE" HMIN,HMAX,EPS,XEND,N,STIFF;
     "BOOLEAN" FIRST,AVAILABLE,STIFF;
     "INTEGER" N;
     "REAL" X,XEND,HMIN,HMAX,EPS;
     "ARRAY" Y,YMAX,SAVE,JACOBIAN;
     "PROCEDURE" DERIV,OUT;
     "CODE" 33080;
 
     MULTISTEP  DELIVERS THE FOLLOWING BOOLEAN VALUE:
             IF DIFFICULTIES  ARE  ENCOUNTERED  DURING  THE  INTEGRATION
             (I.E. SAVE[-1] ^= 0 "OR" SAVE[-2] ^= 0 )  MULTISTEP IS  SET
             TO "FALSE", OTHERWISE MULTISTEP IS SET "TRUE".
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     X:      <VARIABLE>;
             THE INDEPENDENT VARIABLE X.
             CAN BE USED IN DERIV, AVAILABLE ETC.;
             ENTRY: THE INITIAL VALUE X0;
             EXIT : THE FINAL VALUE 'XEND';
     XEND:   <ARITHMETIC EXPRESSION>;
             THE FINAL VALUE OF X ( XEND >= X );
     Y:      <ARRAY IDENTIFIER>;
             "ARRAY" Y[1:6*N];
             THE DEPENDENT VARIABLE;
             ENTRY  Y[1:N] : THE INITIAL VALUES OF THE SOLUTION OF THE
                    SYSTEM OF DIFFERENTIAL EQUATIONS AT X = X0;
             EXIT   Y[1:N] : THE FINAL VALUES OF THE SOLUTION AT X = XEND;
     HMIN,HMAX: <ARITHMETIC EXPRESSION>;
             ENTRY : MINIMAL  RESP. MAXIMAL STEPLENGTH  ALLOWED;
     YMAX:   <ARRAY IDENTIFIER>;
             "ARRAY" YMAX[1:N];
             ENTRY: THE  ABSOLUTE  LOCAL  ERROR  BOUND  DIVIDED  BY  EPS
             EXIT : YMAX[1]  GIVES THE  MAXIMAL VALUE OF THE ENTRY VALUE
                    OF  YMAX[I]  AND  THE  VALUES  OF  ABS(Y[I])  DURING
                    INTEGRATION;
     EPS:    <ARITHMETIC EXPRESSION>;
             THE RELATIVE LOCAL ERROR BOUND;
     FIRST:  <IDENTIFIER>;
             IF FIRST = "TRUE" THEN THE PROCEDURE STARTS THE INTEGRATION
             WITH  A  FIRST  ORDER  ADAMS METHOD AND A STEPLENGTH  EQUAL
             TO HMIN. UPON COMPLETION OF A CALL FIRST:= "FALSE";
             IF   FIRST  =  "FALSE"   THEN   THE   PROCEDURE   CONTINUES
             INTEGRATION;
1SECTION : 5.2.1.1.1.1.F      (DECEMBER 1979)                     PAGE 3
 
 
 
     SAVE:   <ARRAY IDENTIFIER>;
             "ARRAY" SAVE[-38:6*N];
             IN THIS ARRAY THE PROCEDURE STORES INFORMATION WHICH CAN BE
             USED IN A CONTINUING CALL WITH FIRST="FALSE";
             BESIDES THE FOLLOWING MESSAGES ARE DELIVERED:
             SAVE[ 0]=0 : AN ADAMS METHOD HAS BEEN USED;
                      1 : THE PROCEDURE SWITCHED TO GEARS METHOD;
             SAVE[-1]=0 : NO ERROR MESSAGE;
                      1 : WITH  THE HMIN SPECIFIED  THE PROCEDURE CANNOT
                          HANDLE  THE  NONLINEARITY  (DECREASE HMIN!  );
             SAVE[-2]  NUMBER OF TIMES  THAT THE  REQUESTED  LOCAL ERROR
                       BOUND WAS EXCEEDED;
             SAVE[-3]  IF SAVE[-2] IS NONZERO THEN SAVE[-3] GIVES AN
                       ESTIMATE OF THE MAXIMAL LOCAL ERROR BOUND,
                       OTHERWISE SAVE[-3]=0;
     DERIV:  <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:
             "PROCEDURE" DERIV(DF); "ARRAY" DF;
             THIS PROCEDURE SHOULD DELIVER DY[I]/DX IN DF[I];
     AVAILABLE: <BOOLEAN EXPRESSION>;
             IF AN  ANALYTICAL EXPRESSION  OF THE JACOBIAN MATRIX IS NOT
             AVAILABLE THIS EXPRESSION IS SET TO "FALSE";
             OTHERWISE THIS EXPRESSION IS SET TO "TRUE" AND THE
             EVALUATION  OF  THIS  BOOLEAN  EXPRESSION  MUST  EFFECT THE
             FOLLOWING SIDE-EFFECT:
             THE ENTRIES OF THE JACOBIAN MATRIX   D(DY[I]/DX)/DY[J]  ARE
             DELIVERED IN THE ARRAY ELEMENTS  JACOBIAN[I,J];
     JACOBIAN: <ARRAY IDENTIFIER>;
             "ARRAY" JACOBIAN[1:N,1:N];
             AT EACH EVALUATION OF THE BOOLEAN EXPRESSION AVAILABLE WITH
             THE RESULT AVAILABLE:="TRUE", THE JACOBIAN MATRIX HAS TO BE
             ASSIGNED TO THIS ARRAY (SEE THE  EXAMPLE OF USE);
     STIFF:  <BOOLEAN EXPRESSION>;
             IF STIFF = "TRUE"  THE PROCEDURE  SKIPS AN ATTEMPT TO SOLVE
             THE PROBLEM WITH ADAMS-BASHFORTH- OR ADAMS-MOULTON
             METHODS, DIRECTLY USING GEARS METHOD;
     N:      <ARITHMETIC EXPRESSION>;
             THE NUMBER OF EQUATIONS;
     OUT:    <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:
             "PROCEDURE" OUT(H,K);  "VALUE" H,K;  "REAL" H; "INTEGER" K;
             AT THE END OF EACH ACCEPTED STEP OF THE INTEGRATION PROCESS
             THIS PROCEDURE IS CALLED. THE LAST STEPLENGTH USED (H)  AND
             THE ORDER OF THE METHOD (K) ARE DELIVERED.
             AT EACH CALL OF THE  PROCEDURE OUT, THE  CURRENT VALUES  OF
             THE INDEPENDENT VARIABLE (X) AND OF THE SOLUTION (Y[I](X) )
             ARE AVAILABLE  FOR USE. MOREOVER, IN THE  NEIGHBOURHOOD  OF
             THE CURRENT VALUE OF X, ANY VALUE OF  Y[I](X SPECIFIED) CAN
             BE COMPUTED BY MEANS OF THE FOLLOWING INTERPOLATION FORMULA
             Y[I](X SPECIFIED) =
                  SUM(J,0,K, Y[I+J*N] * ((X SPECIFIED - X)/H) ** J ).
 
 
1SECTION : 5.2.1.1.1.1.F      (DECEMBER 1979)                     PAGE 4
 
 
 
 PROCEDURES USED:
 
     MATVEC = CP34011 ,
     DEC    = CP34300 ,
     SOL    = CP34051 .
 
 
 REQUIRED CENTRAL MEMORY: CIRCA  N * ( 2 * N + 5 )  MEMORY PLACES.
 
 
 METHOD AND PERFORMANCE :
     MULTISTEP IS BASED ON TWO LINEAR MULTISTEP METHODS. FOR STIFF
     PROBLEMS IT USES THE BACKWARD DIFFERENTIATION METHODS, FOR
     FOR NON-STIFF PROBLEMS THE ADAMS-BASHFORTH-MOULTON METHODS.
     MULTISTEP IS PROVIDED WITH ORDER, STEPSIZE AND ERROR CONTROL.
 
 
 REFERENCES:
 
     [1].P.W.HEMKER.
         AN ALGOL 60 PROCEDURE  FOR THE SOLUTION  OF STIFF  DIFFERENTIAL
         EQUATIONS.
         MATH. CENTRE, AMSTERDAM. REPORT MR 128/71;
 
 
 EXAMPLE OF USE:
 
     THE SOLUTION AT X=1 AND AT X=10 OF THE DIFFERENTIAL EQUATIONS:
     DY[1]/DX = 0.04 * (1-Y[1]-Y[2]) - Y[1] * ("4*Y[2] + 3"7*Y[1])
     DY[2]/DX = 3"7 Y[1]**2
     WITH THE INITIAL CONDITIONS AT  X = 0 :
     Y[1] = 0  AND Y[2] = 0
     MAY BE OBTAINED BY THE FOLLOWING PROGRAM:
 
1SECTION : 5.2.1.1.1.1.F      (AUGUST 1974)                       PAGE 5
 
 
 
  "BEGIN"
         "BOOLEAN" "PROCEDURE" MULTISTEP(X,XEND,Y,HMIN,HMAX,YMAX,EPS,
                    FIRST,SAVE,DERIV,AVAILABLE,JACOBIAN,STIFF,N,OUT);
         "CODE" 33080;
 
         "BOOLEAN" FIRST;
         "INTEGER" I,J,CF,CJ,CA;
         "REAL" X,XEND,HMIN,EPS,R;
         "ARRAY" Y[1:12],YMAX[1:2],D[-40:12],JAC[1:2,1:2];
 
         "PROCEDURE" DER (F);"ARRAY" F;
         "BEGIN" "REAL" R; CF:=CF+1;
              F[2]:= R:= 3"7*Y[1]*Y[1];
              F[1]:= 0.04*(1-Y[1]-Y[2]) - "4*Y[1]*Y[2] - R;
         "END" F;
 
         "BOOLEAN" "PROCEDURE" AVAIL;
         "BEGIN" "REAL" R; CJ:= CJ+1;
             AVAIL:= "TRUE";
             JAC[2,1]:= R:= 6"7*Y[1];
             JAC[1,1]:= -0.04 - "4*Y[2] - R;
             JAC[1,2]:= -0.04 - "4*Y[1];
             JAC[2,2]:= 0
         "END" JAC AVAIL;
 
         "PROCEDURE" OUT(H,K); "REAL" H; "INTEGER" K; CA:= CA+1;
 
         LABEL:
         OUTPUT(61,"("/,"("HMIN,EPS?")",/")");
         INREAL(60,HMIN); INREAL(60,EPS);
         "IF" HMIN<0 "THEN" "GOTO" ESCAPE;
 
         FIRST:= "TRUE"; CA:=CF:=CJ:=0;
         X:=0; Y[1]:= Y[2]:= 0;
         YMAX[1]:= 0.0001; YMAX[2]:= 1;
 
         "FOR" XEND:= 1, 10 "DO"
         "BEGIN"
             MULTISTEP(X,XEND,Y,HMIN,5,YMAX,EPS,FIRST,D,DER,AVAIL,
                                                JAC,"TRUE",2,OUT);
             OUTPUT(61,"("3(5ZD,2B),2(+.13D"+2D,2B),/")",
                                              CA,CF,CJ,Y[1],Y[2]);
         "END";
 
         "GOTO" LABEL;
         ESCAPE:
     "END"
 
     IT DELIVERS WITH  HMIN = "-10  AND  EPS = "-9:
     240   648   2   +.3074626[602000]"-04   +.3350951[493111]"-01
     315   902   3   +.16233909[62091]"-04   +.15861383[92015]"+00
     (NON-SIGNIFICANT DIGITS ARE PLACED BETWEEN [ ] ).
 
 
1SECTION : 5.2.1.1.1.1.F      (DECEMBER 1979)                     PAGE 6
 
 
 
 SOURCE TEXT(S):
0"CODE"  33080;
  "BOOLEAN" "PROCEDURE" MULTISTEP(X,XEND,Y,HMIN,HMAX,YMAX,EPS,
             FIRST,SAVE,DERIV,AVAILABLE,JACOBIAN,STIFF,N,OUT);
  "VALUE" HMIN,HMAX,EPS,XEND,N,STIFF;
  "BOOLEAN" FIRST,AVAILABLE,STIFF;
  "INTEGER" N;
  "REAL" X,XEND,HMIN,HMAX,EPS;
  "ARRAY" Y,YMAX,SAVE,JACOBIAN;
  "PROCEDURE" DERIV,OUT;
  "BEGIN" "OWN" "BOOLEAN" ADAMS,WITH JACOBIAN;
         "OWN" "INTEGER" M,SAME,KOLD;
         "OWN" "REAL" XOLD,HOLD,A0,TOLUP,TOL,TOLDWN,TOLCONV;
         "BOOLEAN" EVALUATE,EVALUATED,DECOMPOSE,DECOMPOSED,CONV;
         "INTEGER" I,J,L,K,KNEW,FAILS;
         "REAL" H, CH, CHNEW,ERROR,DFI,C;
         "ARRAY" A[0:5],DELTA,LAST DELTA,DF[1:N],JAC[1:N, 1:N],AUX[1:3];
         "INTEGER" "ARRAY" P[1:N];
 
         "REAL" "PROCEDURE" MATVEC(L,U,I,A,B);"CODE" 34011;
         "REAL" "PROCEDURE" DEC(A,N,AUX,P); "CODE" 34300;
         "PROCEDURE" SOL(A,N,P,B); "CODE" 34051;
 
         "REAL" "PROCEDURE" NORM2(AI); "REAL" AI;
         "BEGIN" "REAL" S,A; S:= 1.0"-100;
                 "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" A:= AI/YMAX[I]; S:= S + A * A "END";
                 NORM2:= S
         "END" NORM2;
 
         "PROCEDURE" RESET;
         "BEGIN" "IF" CH < HMIN/HOLD "THEN" CH:= HMIN/HOLD "ELSE"
                 "IF" CH > HMAX/HOLD "THEN" CH:= HMAX/HOLD;
                 X:= XOLD; H:= HOLD * CH; C:= 1;
                 "FOR" J:= 0 "STEP" M "UNTIL" K*M "DO"
                 "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                         Y[J+I]:= SAVE[J+I] * C;
                         C:= C * CH
                 "END";
                 DECOMPOSED:= "FALSE"
         "END" RESET;                                          "COMMENT"
 
1SECTION : 5.2.1.1.1.1.F      (AUGUST 1974)                       PAGE 7
                                                                  ;
 
 
         "PROCEDURE" METHOD;
         "BEGIN" I:= -39;
                 "IF" ADAMS "THEN"
                 "BEGIN" "FOR" C:= 1,1,144,4,0,.5,1,.5,576,144,1,5/12,1,
                               .75,1/6,1436,576,4,.375,1,11/12,1/3,1/24,
                                2844,1436,1,251/720,1,25/24,35/72,
                                5/48,1/120,0,2844,0.1
                         "DO" "BEGIN" I:= I+ 1; SAVE[I]:= C "END"
                 "END" "ELSE"
 
                 "BEGIN" "FOR" C:= 1,1,9,4,0,2/3,1,1/3,36,20.25,1,6/11,
                        1,6/11,1/11,84.028,53.778,0.25,.48,1,.7,.2,.02,
                         156.25, 108.51, .027778, 120/274, 1, 225/274,
                         85/274, 15/274, 1/274, 0, 187.69, .0047361
                         "DO" "BEGIN" I:= I + 1; SAVE[I]:= C "END"
                 "END"
         "END" METHOD;
 
         "PROCEDURE" ORDER;
         "BEGIN" C:= EPS * EPS; J:= (K-1) * (K + 8)/2 - 38;
                 "FOR" I:= 0 "STEP" 1 "UNTIL" K "DO" A[I]:= SAVE[I+J];
                 TOLUP  := C * SAVE[J + K + 1];
                 TOL    := C * SAVE[J + K + 2];
                 TOLDWN := C * SAVE[J + K + 3];
                 TOLCONV:= EPS/(2 * N * (K + 2));
                 A0:= A[0];  DECOMPOSE:= "TRUE";
         "END" ORDER;
 
         "PROCEDURE" EVALUATE JACOBIAN;
         "BEGIN" EVALUATE:= "FALSE";
                 DECOMPOSE:= EVALUATED:= "TRUE";
                 "IF" AVAILABLE "THEN" "ELSE"
                 "BEGIN" "REAL" D; "ARRAY" FIXY,FIXDY,DY[1:N];
                         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                         FIXY[I]:= Y[I];
                         DERIV(FIXDY);
                         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
                         "BEGIN" D:= "IF" EPS > ABS(FIXY[J])
                                 "THEN" EPS * EPS
                                 "ELSE" EPS * ABS(FIXY[J]);
                                 Y[J]:= Y[J] + D; DERIV(DY);
                                 "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                                 JACOBIAN[I,J]:= (DY[I]-FIXDY[I])/D;
                                 Y[J]:= FIXY[J]
                         "END"
                 "END"
         "END" EVALUATE JACOBIAN;                              "COMMENT"
 
1SECTION : 5.2.1.1.1.1.F      (DECEMBER 1979)                     PAGE 8
                                                                  ;
 
 
         "PROCEDURE" DECOMPOSE JACOBIAN;
         "BEGIN" DECOMPOSE:= "FALSE";
                 DECOMPOSED:= "TRUE"; C:= -A0 * H;
                 "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                         JAC[I,J]:= JACOBIAN[I,J] * C;
                         JAC[J,J]:= JAC[J,J] + 1
                 "END";
                 AUX[2]:=1.0"-12;
                 DEC(JAC,N,AUX,P)
         "END" DECOMPOSE JACOBIAN;
 
         "PROCEDURE" CALCULATE STEP AND ORDER;
         "BEGIN" "REAL" A1,A2,A3;
                 A1:= "IF" K <= 1 "THEN" 0 "ELSE"
                      0.75 * (TOLDWN/NORM2(Y[K*M+I])) ** (0.5/K);
                 A2:= 0.80 * (TOL/ERROR) ** (0.5/(K + 1));
                 A3:= "IF" K >= 5 "OR" FAILS ^= 0
                      "THEN" 0 "ELSE"
                      0.70 * (TOLUP/NORM2(DELTA[I] - LAST DELTA[I])) **
                      (0.5/(K+2));
 
                 "IF" A1 > A2 "AND" A1 > A3 "THEN"
                 "BEGIN" KNEW:= K-1; CHNEW:= A1 "END" "ELSE"
                 "IF" A2 > A3 "THEN"
                 "BEGIN" KNEW:= K  ; CHNEW:= A2 "END" "ELSE"
                 "BEGIN" KNEW:= K+1; CHNEW:= A3 "END"
         "END" CALCULATE STEP AND ORDER;
 
         "IF" FIRST "THEN"
         "BEGIN" FIRST:= "FALSE"; M:= N;
                 "FOR" I:= -1,-2,-3 "DO" SAVE[I]:= 0;
                 OUT(0,0);
                 ADAMS:= "NOT" STIFF; WITH JACOBIAN:= "NOT" ADAMS;
                 "IF" WITH JACOBIAN "THEN" EVALUATE JACOBIAN;
                 METHOD;
         NEW START: K:= 1; SAME:= 2; ORDER; DERIV(DF);
                 H:= "IF" "NOT" WITH JACOBIAN "THEN" HMIN "ELSE"
                 SQRT(2 * EPS/SQRT(NORM2 (MATVEC(1,N,I,JACOBIAN,DF))));
                 "IF" H > HMAX "THEN" H:= HMAX "ELSE"
                 "IF" H < HMIN "THEN" H:= HMIN;
                 XOLD:= X; HOLD:= H; KOLD:= K; CH:= 1;
                 "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" SAVE[I]:= Y[I]; SAVE[M+I]:= Y[M+I]:= DF[I] * H
                 "END";
                 OUT(0,0)
         "END" "ELSE"
         "BEGIN" WITH JACOBIAN:= "NOT" ADAMS; CH:= 1;
                 K:=KOLD; RESET; ORDER;
                 DECOMPOSE:= WITH JACOBIAN
         "END";
         FAILS:= 0;                                            "COMMENT"
 
1SECTION : 5.2.1.1.1.1.F      (AUGUST 1974)                       PAGE 9
                                                                  ;
 
 
         "FOR" L:= 0 "WHILE" X < XEND "DO"
         "BEGIN" "IF" X + H <= XEND "THEN" X:= X + H "ELSE"
                 "BEGIN" H:= XEND-X; X:= XEND; CH:= H/HOLD; C:= 1;
                         "FOR" J:= M "STEP" M "UNTIL" K*M "DO"
                         "BEGIN" C:= C* CH;
                                 "FOR" I:= J+1 "STEP" 1 "UNTIL" J+N "DO"
                                 Y[I]:= Y[I] * C
                         "END";
                         SAME:= "IF" SAME<3 "THEN" 3 "ELSE" SAME+1;
                 "END";
 
                 "COMMENT" PREDICTION;
                 "FOR" L:= 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" "FOR" I:= L "STEP" M "UNTIL" (K-1)*M+L "DO"
                         "FOR" J:= (K-1)*M+L "STEP" -M "UNTIL" I "DO"
                         Y[J]:= Y[J] + Y[J+M];
                         DELTA[L]:= 0
                 "END";  EVALUATED:= "FALSE";
 
         "COMMENT" CORRECTION AND ESTIMATION LOCAL ERROR;
         "FOR" L:= 1,2,3 "DO"
         "BEGIN" DERIV(DF);
                 "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                 DF[I]:= DF[I] * H - Y[M+I];
                 "IF" WITH JACOBIAN "THEN"
                 "BEGIN" "IF" EVALUATE "THEN" EVALUATE JACOBIAN;
                         "IF" DECOMPOSE "THEN" DECOMPOSE JACOBIAN;
                         SOL(JAC,N,P,DF)
                 "END";
 
                 CONV:= "TRUE";
                 "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" DFI:= DF[I];
                         Y[  I]:= Y[  I] + A0 * DFI;
                         Y[M+I]:= Y[M+I] +      DFI;
                         DELTA[I]:= DELTA[I] +  DFI;
                         CONV:= CONV "AND" ABS(DFI) < TOLCONV * YMAX[I]
                 "END";
                 "IF" CONV "THEN"
                 "BEGIN" ERROR:= NORM2(DELTA[I]);
                         "GOTO" CONVERGENCE
                 "END"
         "END";                                                "COMMENT"
 
1SECTION : 5.2.1.1.1.1.F      (DECEMBER 1979)                    PAGE 10
                                                                  ;
 
 
         "COMMENT" ACCEPTANCE OR REJECTION;
         "IF" "NOT" CONV "THEN"
         "BEGIN" "IF" "NOT" WITH JACOBIAN "THEN"
                 "BEGIN" EVALUATE:= WITH JACOBIAN:= SAME >= K
                            "OR" H<1.1 * HMIN;
                         "IF" "NOT" WITH JACOBIAN "THEN" CH:= CH/4;
                 "END" "ELSE"
                 "IF" "NOT" DECOMPOSED "THEN" DECOMPOSE:= "TRUE" "ELSE"
                 "IF" "NOT" EVALUATED  "THEN" EVALUATE := "TRUE" "ELSE"
                 "IF" H > 1.1 * HMIN   "THEN" CH:= CH/4 "ELSE"
                 "IF" ADAMS            "THEN" "GOTO" TRY CURTISS "ELSE"
                 "BEGIN" SAVE[-1]:= 1; "GOTO" RETURN "END";
 
                 RESET
         "END" "ELSE" CONVERGENCE:
 
         "IF" ERROR > TOL "THEN"
         "BEGIN" FAILS:= FAILS + 1;
                 "IF" H > 1.1 * HMIN "THEN"
                 "BEGIN" "IF" FAILS > 2 "THEN"
                         "BEGIN" "IF" ADAMS "THEN"
                                 "BEGIN" ADAMS:= "FALSE"; METHOD "END";
                                 KOLD:= 0; RESET; "GOTO" NEW START
                         "END" "ELSE"
                         "BEGIN" CALCULATE STEP AND ORDER;
                                 "IF" KNEW ^= K "THEN"
                                 "BEGIN" K:= KNEW; ORDER "END";
                                 CH:= CH * CHNEW; RESET
                         "END"
                 "END" "ELSE"
                 "BEGIN" "IF" ADAMS "THEN" TRY CURTISS:
                         "BEGIN" ADAMS:= "FALSE"; METHOD
                         "END" "ELSE"
                         "IF" K = 1 "THEN"
                         "BEGIN" "COMMENT" VIOLATE EPS CRITERION;
                                 C:= EPS * SQRT(ERROR/TOL);
                                 "IF" C > SAVE[-3] "THEN" SAVE[-3]:= C;
                                 SAVE[-2]:= SAVE[-2] + 1;
                                 SAME:= 4; "GOTO" ERROR TEST OK
                         "END";
                         K:= KOLD:= 1; RESET; ORDER; SAME:= 2
                 "END"
         "END" "ELSE" ERROR TEST OK:
         "BEGIN"                                               "COMMENT"
1SECTION : 5.2.1.1.1.1.F      (AUGUST 1974)                      PAGE 11
                                                                  ;
 
 
                 FAILS:= 0;
                 "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" C:= DELTA[I];
                         "FOR" L:= 2 "STEP" 1 "UNTIL" K "DO"
                         Y[L*M+I]:= Y[L*M+I] + A[L] * C;
                         "IF" ABS(Y[I]) > YMAX[I] "THEN"
                              YMAX[I]:=  ABS(Y[I])
                 "END";
 
                 SAME:= SAME-1;
                 "IF" SAME= 1 "THEN"
                 "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
                         LAST DELTA[I]:= DELTA[I]
                 "END" "ELSE"
                 "IF" SAME= 0 "THEN"
                 "BEGIN" CALCULATE STEP AND ORDER;
                         "IF" CHNEW > 1.1 "THEN"
                         "BEGIN" DECOMPOSED:= "FALSE";
                                 "IF" K ^= KNEW "THEN"
                                 "BEGIN" "IF" KNEW > K "THEN"
                                         "BEGIN" "FOR" I:= 1 "STEP" 1
                                              "UNTIL" N "DO" Y[KNEW*M+I]
                                                 := DELTA[I] * A[K]/KNEW
                                         "END";
                                         K:= KNEW; ORDER
                                 "END";
                                 SAME:= K+1;
                                 "IF" CHNEW * H > HMAX
                                     "THEN" CHNEW:= HMAX/H;
                                 H:= H * CHNEW; C:= 1;
                                 "FOR" J:= M "STEP" M "UNTIL" K*M "DO"
                                 "BEGIN" C:= C * CHNEW;
                                         "FOR" I:= J+1 "STEP" 1 "UNTIL"
                                         J+N "DO" Y[I]:= Y[I] * C
                                 "END"
                         "END"
                         "ELSE" SAME:= 10
                 "END";
                 "IF" X ^= XEND "THEN"
                 "BEGIN" XOLD:= X; HOLD:= H; KOLD:= K; CH:= 1;
                         "FOR" I:= K * M + N "STEP" -1 "UNTIL" 1 "DO"
                         SAVE[I]:= Y[I];
                         OUT(H,K)
                 "END"
         "END" CORRECTION AND ESTIMATION LOCAL ERROR;
         "END" STEP;
 
         RETURN: SAVE[0]:= "IF" ADAMS "THEN" 0 "ELSE" 1;
         MULTISTEP:= SAVE[-1]= 0 "AND" SAVE[-2]=0
  "END" MULTISTEP;
         "EOP"
1SECTION : 5.2.1.1.1.1.G      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : DIFFSYS.       
        
        
 AUTHORS   : R.BULIRSCH AND J.STOER.  
        
        
 CONTRIBUTOR: K.DEKKER.     
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731231.
        
        
 BRIEF DESCRIPTION: 
        
     DIFFSYS SOLVES AN INITIAL VALUE PROBLEM ,FOR A SYSTEM OF FIRST 
     ORDER ORDINARY DIFFERENTIAL EQUATIONS DY / DX = F(X,Y).
     THE METHOD IS RECOMMENDED IF HIGH ACCURACY IS DESIRED. 
     DIFFSYS IS NOT SUITED FOR STIFF EQUATIONS. 
        
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEMS,
     SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS. 
        
        
1SECTION : 5.2.1.1.1.1.G      (FEBRUARY 1979)                     PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE DIFFSYS READS:
     "PROCEDURE" DIFFSYS(X,XE,N,Y,DERIVATIVE,AETA,RETA,S,H0,OUTPUT);
     "VALUE" N;   
     "INTEGER" N; 
     "REAL" X,XE,AETA,RETA,H0;
     "ARRAY" Y,S; 
     "PROCEDURE" DERIVATIVE,OUTPUT;   
     "CODE" 33180;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:      <VARIABLE>;    
             THE INDEPENDENT VARIABLE  ;
             ENTRY: THE INITIAL VALUE X0;       
             EXIT : THE FINAL VALUE XE; 
     XE:     <ARITHMETIC EXPRESSION>; 
             THE FINAL VALUE OF X (XE>=X);      
     N:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF EQUATIONS; 
     Y:      <ARRAY IDENTIFIER>;      
             "REAL" "ARRAY" Y[1:N];   
             THE DEPENDENT VARIABLE;  
             ENTRY: THE  INITIAL  VALUES  OF  THE SYSTEM OF DIFFERENTIAL      
                    EQUATIONS AT X=X0;
             EXIT : THE FINAL VALUES OF THE SOLUTION AT X=XE;       
     DERIVATIVE: <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS: 
             "PROCEDURE" DERIVATIVE(X,Y,DY); "REAL" X; "ARRAY" Y,DY;
             THIS  PROCEDURE  SHOULD  DELIVER  THE  RIGHT  HAND  SIDE OF      
             THE I-TH DIFFERENTIAL EQUATION AT THE POINT (X,Y) AS DY[I],      
             I=1,...,N;     
     AETA:   <ARITHMETIC EXPRESSION>; 
             REQUIRED  ABSOLUTE  PRECISION  IN THE  INTEGRATION PROCESS;      
     RETA:   <ARITHMETIC EXPRESSION>; 
             REQUIRED  RELATIVE  PRECISION  IN THE  INTEGRATION PROCESS;      
     S:      <ARRAY IDENTIFIER>;      
             "REAL" "ARRAY" S[1:N];   
             THE ARRAY S IS USED TO CONTROL THE ACCURACY OF THE COMPUTED      
             VALUES OF Y;   
             ENTRY: IT IS ADVISABLE TO SET S[I]=0, I=1,...,N;       
             EXIT : THE MAXIMUM VALUE  OF ABS(Y[I]), ENCOUNTERED  DURING      
                    INTEGRATION, IF THIS VALUE EXCEEDS THE VALUE OF S[I]      
                    ON ENTRY; 
     H0:     <VARIABLE>;    
             THE INITIAL STEP TO BE TAKEN;      
     OUTPUT: <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS; 
             "PROCEDURE" OUTPUT;      
             THIS  PROCEDURE  IS  CALLED  AT THE END OF EACH INTEGRATION      
             STEP ; THE USER CAN ASK FOR OUTPUT OF SOME PARAMETERS , FOR      
             EXAMPLE X, Y, S. 
        
        
1SECTION : 5.2.1.1.1.1.G      (FEBRUARY 1979)                     PAGE 3      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 REQUIRED CENTRAL MEMORY : CIRCA 28* N  MEMORY PLACES.    
        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE  DIFFSYS  IS  A SLIGHT  MODIFICATION OF THE ALGORITHM      
     PUBLISHED BY BULIRSCH AND STOER (SEE REF[1]) . BY THIS MODIFICATION      
     INTEGRATION  FROM X0  UNTIL XE  CAN BE  PERFORMED  BY  ONE  CALL OF      
     DIFFSYS. A NUMBER OF INTEGRATION STEPS ARE TAKEN, STARTING WITH THE      
     INITIAL STEP H0. IN EACH INTEGRATION STEP A NUMBER OF SOLUTIONS ARE      
     COMPUTED BY MEANS OF THE MODIFIED MIDPOINT RULE . EXTRAPOLATION  IS      
     USED TO IMPROVE THESE SOLUTIONS,UNTIL THE REQUIRED ACCURACY IS MET.      
     AN INTEGRATION STEP  IS REJECTED, IF THE ACCURACY REQUIREMENTS  ARE      
     NOT FULFILLED  AFTER NINE EXTRAPOLATION STEPS . IN THESE CASES  THE      
     INTEGRATION STEP IS REJECTED , AND INTEGRATION IS TRIED  AGAIN WITH      
     THE INTEGRATION STEP HALVED.     
     THE ALGORITHM IS FOR EACH STEP A VARIABLE ORDER METHOD (THE HIGHEST      
     ORDER IS 14 ), AND USES A VARIABLE NUMBER OF  FUNCTION EVALUATIONS,      
     DEPENDING ON THE ORDER (MINIMUM IS 3, MAXIMUM IS 217). 
     THE ALGORITHM IS LESS SENSITIVE TO TOO SMALL VALUES  OF THE INITIAL      
     STEPSIZE THAN THE ORIGINAL ALGORITHM ( SEE REF [2] );  HOWEVER  BAD      
     GUESSES REQUIRE STILL SOME MORE COMPUTATIONS.
        
 REFERENCES:      
        
     [1]. R.BULIRSCH AND J.STOER.     
          NUMERICAL  TREATMENT  OF  ORDINARY  DIFFERENTIAL  EQUATIONS BY      
          EXTRAPOLATION METHODS.      
          NUMERISCHE MATHEMATIK, VOLUME 8, PAGE 1-13, 1965. 
        
     [2]. PHYLLIS FOX.      
          A  COMPARATIVE  STUDY  OF  COMPUTER  PROGRAMS  FOR INTEGRATING      
          DIFFERENTIAL EQUATIONS.     
          COMMUNICATIONS OF THE A.C.M., VOLUME 15, PAGE 941-948, 1972.
        
     [3]. T.E.HULL, W.H.ENRIGHT, B.M.FELLEN AND A.E.SEDGWICK.       
          COMPARING   NUMERICAL   METHODS  FOR   ORDINARY   DIFFERENTIAL      
          EQUATIONS.
          SIAM JOURNAL ON NUMERICAL ANALYSIS,VOLUME 9,PAGE 603-635,1972.      
        
1SECTION : 5.2.1.1.1.1.G      (AUGUST 1974)                       PAGE 4      
        
        
        
 EXAMPLE OF USE:  
        
     THE FOLLOWING  PROGRAM  ILLUSTRATES  THE COSTS  AND  THE ACCURACIES      
     WHICH ARE OBTAINED WHEN SOLVING A SYSTEM OF DIFFERENTIAL  EQUATIONS      
     ARISING FROM THE RESTRICTED PROBLEM OF THREE BODIES ( SEE REF[1] ).      
     THE SOLUTION IS A CLOSED ORBIT WITH PERIOD T=6.192169331396.   
        
     "BEGIN"      
         "PROCEDURE" DIFFSYS(X,XE,N,Y,DERIVATIVE,AETA,RETA,S,H0,OUTPUT);      
         "CODE" 33180;      
         "INTEGER" PASSES,K;
         "REAL" X,XE,TIME,TOL,H0;     
         "REAL" "ARRAY" Y,S[1:4];     
        
         "PROCEDURE" DER(X,Y,DY); "REAL" X; "ARRAY" Y,DY; 
         "BEGIN" "REAL" MU,MU1,Y1,Y2,Y3,Y4,S1,S2; 
             MU:=1/82.45; MU1:=1-MU;  
             PASSES:=PASSES+1;
             Y1:=Y[1]; Y2:=DY[1]:=Y[2]; Y3:=Y[3]; Y4:=DY[3]:=Y[4];  
             S1:=(Y1+MU)**2+Y3**2; S2:=(Y1-MU1)**2+Y3**2; 
             S1:=S1*SQRT(S1); S2:=S2*SQRT(S2);  
             DY[2]:=Y1+2*Y4-MU1*(Y1+MU)/S1-MU*(Y1-MU1)/S2;
             DY[4]:=Y3-2*Y2-MU1*Y3/S1-MU*Y3/S2  
         "END";   
        
         "PROCEDURE" OUT;   
         "BEGIN" K:=K+1;    
             "IF" X>=XE "THEN"
             OUTPUT(61,"("2(-5ZD),2(4B+Z.3DB3DB3DB3D),-5ZD.3D,/")",K, 
             PASSES,Y[1],Y[3],CLOCK-TIME)       
         "END";   
        
         OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
         "("PRINTED BY THIS PROGRAM")",//,      
         "("    THE RESULTS WITH DIFFSYS - H0=.2 - ARE: ")",/,      
         "("     K   DER.EV.         Y[1]                 Y[3]      ")",      
         "("       TIME")",/")");     
         "FOR" TOL:="-4,"-6,"-8,"-10,"-12 "DO"  
         "BEGIN" PASSES:=K:=0; X:=0; XE:=6.192169331396;  
             Y[1]:=1.2; Y[2]:=Y[3]:=0; Y[4]:=-1.04935750983;
             S[1]:=S[2]:=S[3]:=S[4]:=0; H0:=.2; TIME:=CLOCK;
             DIFFSYS(X,XE,4,Y,DER,TOL,TOL,S,H0,OUT);      
         "END"    
     "END"
        
     THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:   
        
     THE RESULTS WITH DIFFSYS - H0=.2 - ARE:    
      K   DER.EV.         Y[1]                  Y[3]             TIME 
      30   2591    +1.320 357 347 741     -.032 645 454 836      5.686
      33   3414    +1.200 078 037 878     -.000 053 906 067      7.455
      37   4213    +1.200 003 282 801     -.000 002 363 741      9.267
      44   4618    +1.199 999 999 711     -.000 000 000 095     10.242
      56   6299    +1.200 000 000 003     -.000 000 000 090     13.827
        
        
1SECTION : 5.2.1.1.1.1.G      (AUGUST 1974)                       PAGE 5      
        
        
        
 SOURCE TEXT:     
0"CODE"" 33180;   
 "PROCEDURE" DIFFSYS(X,XE,N,Y,DERIVATIVE,AETA,RETA,S,H0,OUTPUT);    
 "VALUE" N;       
 "INTEGER" N;     
 "REAL" X,XE,AETA,RETA,H0;  
 "ARRAY" Y,S;     
 "PROCEDURE" DERIVATIVE,OUTPUT;       
 "BEGIN" "REAL" A,B,B1,C,G,H,U,V,TA,FC; "INTEGER" I,J,K,KK,JJ,L,M,R,SR;       
     "ARRAY" YA,YL,YM,DY,DZ[1:N],DT[1:N,0:6],D[0:6],YG,YH[0:7,1:N]; 
     "BOOLEAN" KONV,B0,BH,LAST;       
     LAST:="FALSE"; H:=H0;  
 NEXT: "IF" H*1.1>=XE-X "THEN"
     "BEGIN" LAST:="TRUE"; H0:=H; H:=XE-X+"-13 "END";     
     DERIVATIVE(X,Y,DZ); BH:="FALSE"; 
     "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" YA[I]:=Y[I];      
 ANF: A:=H+X; FC:=1.5; B0:="FALSE"; M:=1; R:=2; SR:=3; JJ:=-1;      
     "FOR" J:=0 "STEP" 1 "UNTIL" 9 "DO" 
     "BEGIN" "IF" B0 "THEN" 
         "BEGIN" D[1]:=16/9; D[3]:=64/9; D[5]:=256/9 "END"
         "ELSE" "BEGIN" D[1]:=9/4; D[3]:=9; D[5]:=36 "END"; 
         KONV:="TRUE";      
         "IF" J>6 "THEN" "BEGIN" L:=6; D[6]:=64; FC:=.6*FC "END"    
         "ELSE" "BEGIN" L:=J; D[L]:=M*M "END";  
         M:=M*2; G:=H/M; B:=G*2;      
         "IF" BH "AND" J<8 "THEN"     
         "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"       
             "BEGIN" YM[I]:=YH[J,I]; YL[I]:=YG[J,I] "END" 
         "END"    
         "ELSE"   
         "BEGIN"                                               "COMMENT"      
1SECTION : 5.2.1.1.1.1.G      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
                 KK:=(M-2)/2; M:=M-1; 
             "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" 
             "BEGIN" YL[I]:=YA[I]; YM[I]:=YA[I]+G*DZ[I] "END";      
             "FOR" K:=1 "STEP" 1 "UNTIL" M "DO" 
             "BEGIN" DERIVATIVE(X+K*G,YM,DY);   
                 "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"       
                 "BEGIN" U:=YL[I]+B*DY[I]; YL[I]:=YM[I]; YM[I]:=U;  
                     U:=ABS(U); "IF" U>S[I] "THEN" S[I]:=U
                 "END";     
                 "IF" K=KK "AND" K^=2 "THEN"    
                 "BEGIN" JJ:=JJ+1; "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" 
                     "BEGIN" YH[JJ,I]:=YM[I]; YG[JJ,I]:=YL[I] "END" 
                 "END"      
             "END"
         "END";   
         DERIVATIVE(A,YM,DY); 
         "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"     
         "BEGIN" V:=DT[I,0]; TA:=C:=DT[I,0]:=(YM[I]+YL[I]+G*DY[I])/2; 
             "FOR" K:=1 "STEP" 1 "UNTIL" L "DO" 
             "BEGIN" B1:=D[K]*V; B:=B1-C; U:=V; 
                 "IF" B^=0 "THEN"     
                 "BEGIN" B:=(C-V)/B; U:=C*B; C:=B1*B "END"; 
                 V:=DT[I,K]; DT[I,K]:=U; TA:=U+TA 
             "END"; 
             "IF" ABS(Y[I]-TA)>RETA*S[I]+AETA "THEN" KONV:="FALSE"; 
             Y[I]:=TA       
         "END";   
         "IF" KONV "THEN" "GOTO" END; 
         D[2]:=4; D[4]:=16; B0:=^B0; M:=R; R:=SR; SR:=M*2 
     "END";       
     BH:=^BH; LAST:="FALSE"; H:=H/2; "GOTO" ANF;
 END: H:=FC*H; X:=A; OUTPUT; "IF" "NOT" LAST "THEN" "GOTO" NEXT;    
 "END" DIFFSYS;   
         "EOP"    
1SECTION : 5.2.1.1.1.1.H     (FEBRUARY 1979)                     PAGE 1
 
 
 
 PROCEDURE : ARK.
 
 
 AUTHOR:  P.A. BEENTJES.
 
 
 INSTITUTE:  MATHEMATICAL CENTRE.
 
 
 RECEIVED: 740510.
 
 
 BRIEF DESCRIPTION:
 
     ARK SOLVES AN INITIAL VALUE PROBLEM, FOR A SYSTEM  OF FIRST ORDER
     ORDINARY DIFFERENTIAL EQUATIONS . ARK IS RECOMMENDED FOR THE
     INTEGRATION OF SEMI-DISCRETE PARABOLIC AND HYPERBOLIC
     INITIAL-BOUNDARY PROBLEMS.
 
 
 KEYWORDS:
     INITIAL VALUE PROBLEM,
     SEMI-DISCRETE PARABOLIC AND HYPERBOLIC PROBLEM.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" ARK (T, TE, M0, M, U, DERIVATIVE, DATA, OUT);
     "INTEGER" M0, M; "REAL" T, TE; "ARRAY" U, DATA;
     "PROCEDURE" DERIVATIVE, OUT;
     "CODE" 33061;
 
     ARK :   INTEGRATES THE SYSTEM OF ORDINARY DIFFERENTIAL  EQUATIONS
             DU / DT = H(T, U),  U = U0   AT   T = T0.
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     T:      <VARIABLE>;
             THE  INDEPENDENT  VARIABLE  T;
             ENTRY: THE INITIAL VALUE T0;
             EXIT : THE FINAL VALUE TE;
     TE:     <ARITHMETIC EXPRESSION>;
             THE FINAL VALUE OF T (TE >= T);
     M0,M:   <ARITHMETIC EXPRESSION>;
             INDICES  OF  THE  FIRST  AND  LAST EQUATION OF THE SYSTEM;
     U:      <ARRAY IDENTIFIER>;
             "ARRAY" U[M0 : M];
             ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM OF
                    DIFFERENTIAL EQUATIONS AT T = T0;
             EXIT : THE VALUES OF THE SOLUTION AT T = TE;
1SECTION : 5.2.1.1.1.1.H      (FEBRUARY 1979)                     PAGE 2
 
 
 
     DERIVATIVE: <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:
             "PROCEDURE" DERIVATIVE(T, V); "REAL" T; "ARRAY" V;
             THIS PROCEDURE  PERFORMS  AN  EVALUATION OF THE RIGHT HAND
             SIDE  OF THE SYSTEM WITH DEPENDENT VARIABLES V[M0 : M] AND
             INDEPENDENT VARIABLE T; UPON COMPLETION OF DERIVATIVE, THE
             RIGHT HAND SIDE   SHOULD  BE   OVERWRITTEN  ON  V[M0 : M];
     DATA:   <ARRAY IDENTIFIER>;
             "ARRAY" DATA[1 : 10 + DATA[1]];
             IN ARRAY DATA ONE SHOULD GIVE:
             DATA[1]: THE   NUMBER   OF   EVALUATIONS  OF   H(T, U) PER
                      INTEGRATION STEP(DATA[1] >= DATA[2]);
             DATA[2]: THE ORDER OF ACCURACY OF THE METHOD (DATA[2]<=3);
             DATA[3]: STABILITY BOUND(SEE REFERENCE [3]);
             DATA[4]: THE SPECTRAL RADIUS  OF THE JACOBIAN MATRIX  WITH
                      RESPECT TO  THOSE EIGENVALUES, WHICH ARE  LOCATED
                      IN THE NON-POSITIVE HALF PLANE;
             DATA[5]: THE MINIMAL STEPSIZE;
             DATA[6]: THE ABSOLUTE TOLERANCE;
             DATA[7]: THE RELATIVE TOLERANCE;
                      IF  BOTH  DATA[6]  AND  DATA[7] ARE NEGATIVE, THE
                      INTEGRATION  IS  PERFORMED  WITH  A CONSTANT STEP
                      DATA[5];
             DATA[8]: DATA[8]  SHOULD  BE  0  IF  ARK  IS  CALLED   FOR
                      A FIRST TIME;  FOR  CONTINUED  INTEGRATION
                      DATA[8] SHOULD NOT BE CHANGED;
             DATA[11], ..., DATA[10 + DATA[1]]: POLYNOMIAL COEFFICIENTS
                      (SEE REFERENCE [3]);
             AFTER  EACH  STEP THE FOLLOWING BY-PRODUCTS ARE DELIVERED:
             DATA[8]: THE   NUMBER  OF   INTEGRATION  STEPS  PERFORMED;
             DATA[9]: AN  ESTIMATION  OF THE  LOCAL  ERROR  LAST  MADE;
             DATA[10]: INFORMATIVE MESSAGES:
                      DATA[10] = 0: NO DIFFICULTIES;
                      DATA[10] = 1: MINIMAL  STEPLENGTH   EXCEEDS   THE
                              STEPLENGTH   PRESCRIBED   BY    STABILITY
                              THEORY, I.E. DATA[5] > DATA[3] / DATA[4];
                              (TERMINATION OF ARK);
                              DECREASE  MINIMAL  STEPLENGTH;
             IF NECESSARY, DATA[I],I = 4(1)7, CAN BE UPDATED(AFTER EACH
             STEP) BY MEANS OF PROCEDURE OUT;
     OUT:    <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:
             "PROCEDURE" OUT;
             AFTER  EACH  INTEGRATION STEP PERFORMED INFORMATION CAN BE
             OBTAINED  OR UPDATED BY THIS PROCEDURE, E.G. THE VALUES OF
             T, U[M0 : M] AND DATA[I], I = 4(1)10.
1SECTION : 5.2.1.1.1.1.H      (FEBRUARY 1979)                     PAGE 3
 
 
 
 DATA AND RESULTS:
 
     FOR THE INDICES  M0  AND  M  THE  FOLLOWING  REMARKS  CAN BE MADE:
     WHEN THE METHOD OF LINES  IS  APPLIED  TO  HYPERBOLIC DIFFERENTIAL
     EQUATIONS THE NUMBER OF RELEVANT ORDINARY  DIFFERENTIAL  EQUATIONS
     DECREASES  DURING  THE  INTEGRATION  PROCESS;  IN   PROCEDURE  ARK
     THIS   MAY   BE   REALIZED   BY INTEGERS  M0  AND  M,  WHICH   ARE
     DEFINED AS FUNCTIONS OF THE NUMBER OF RIGHT HAND SIDE EVALUATIONS.
     A SELECTION OF POSSIBLE ENTRIES FOR ARRAY DATA (DEPENDENT ON THE
     KIND OF INITIAL VALUE PROBLEM) IS GIVEN IN REFERENCE [4],SECTION 8.
 
 
 PROCEDURES USED:
 
     INIVEC = CP31010,
     MULVEC = CP31020,
     DUPVEC = CP31030,
     VECVEC = CP34010,
     ELMVEC = CP34020,
     DECSOL = CP34301.
 
 
 REQUIRED CENTRAL MEMORY : CIRCA 75 + 2 * (M - M0) MEMORY PLACES.
 
 
 METHOD AND PERFORMANCE:
 
     ARK IS  AN  IMPLEMENTATION  OF  LOW  ORDER  STABILIZED RUNGE KUTTA
     METHODS (SEE REFERENCE [1]);
     AUTOMATIC STEPSIZE CONTROL IS PROVIDED BUT STEP-REJECTION HAS BEEN
     EXCLUDED IN ORDER TO SAVE STORAGE;
     BECAUSE OF ITS LIMITED STORAGE REQUIREMENTS AND ADAPTIVE STABILITY
     FACILITIES THE METHOD IS WELL SUITED FOR THE SOLUTION  OF  INITIAL
     BOUNDARY  VALUE  PROBLEMS   FOR  PARTIAL  DIFFERENTIAL  EQUATIONS;
     NUMERICAL   RESULTS,   OBTAINED   WITH   A    SLIGHTLY   DIFFERENT
     IMPLEMENTATION CAN BE FOUND IN REFERENCE [2].
 
 
 REFERENCES:
 
     [1]. P.J. VAN DER HOUWEN.
             STABILIZED RUNGE KUTTA METHOD WITH LIMITED
             STORAGE REQUIREMENTS.
             MATH. CENTR. REPORT  TW 124/71;
 
     [2]. P.A. BEENTJES.
             AN ALGOL 60 VERSION OF STABILIZED RUNGE KUTTA
             METHODS (DUTCH).
             MATH. CENTR. REPORT NR 23/72;
1SECTION : 5.2.1.1.1.1.H      (DECEMBER 1975)                     PAGE 4
 
 
 
     [3]. P.J. VAN DER HOUWEN, J. KOK.
             NUMERICAL SOLUTION OF A MINIMAX PROBLEM.
             MATH. CENTR. REPORT TW 123/71;
 
     [4]. P.J. VAN DER HOUWEN ET AL.
             ONE STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS, I.I.I.
             NUMERICAL EXAMPLES,
             MATH. CENTR. REPORT TW 130/71.
 
 
 EXAMPLE OF USE:
 
     THE VALUES OF
 
     1. Y(1) AND Y(2) OF THE INITIAL VALUE PROBLEM
        DY / DX = Y - 2 * X / Y,    Y(0) = 1
 
      AND
 
     2. U(.6, 0) OF THE CAUCHY PROBLEM (SEE REFERENCE [2]):
        DU / DT = .5 * DU / DX,  U(0, X) = EXP(-X * X)
 
     MAY BE OBTAINED BY THE FOLLOWING PROGRAM:
 
 
 "BEGIN" "INTEGER" M0, M, I; "REAL" T, TE, DAT;
     "ARRAY" Y[1 : 1], U[-150 : 150], DATA[1 : 14];
 
     "PROCEDURE" ARK
     (T, TE, M0, M, U, DERIVATIVE, DATA, OUT); "CODE" 33061;
 
     "PROCEDURE" DER1(T, V); "REAL" T; "ARRAY" V;
     V[1]:= V[1] - 2 * T / V[1];
 
     "PROCEDURE" DER2(T, V); "REAL" T; "ARRAY" V;
     "BEGIN" "INTEGER" J; "REAL" V1, V2, V3;
        V2:= V[M0]; M0:= M0 + 1; M:= M - 1; V3:= V[M0];
        "FOR" J:= M0 "STEP" 1 "UNTIL" M "DO"
        "BEGIN" V1:= V2; V2:= V3; V3:= V[J + 1];
            V[J]:= 250 * (V3 - V1) / 3
        "END"
     "END" DER2;
 
1SECTION : 5.2.1.1.1.1.H     (DECEMBER 1975)                     PAGE 5
 
 
 
     "PROCEDURE" OUT1;
     "IF" T = TE "THEN"
     "BEGIN" "IF" T = 1 "THEN" OUTPUT(61, "("/, "(" PROBLEM 1")", //,
         "(" X  NUMBER OF INTEGRATION STEPS Y(COMPUTED)  Y(EXACT)")",
         //")");
         OUTPUT(61, "("ZD, 13ZD,12B,2(-3ZD.7D),"("...")", /")",
         T, DATA[8], Y[1], SQRT(2 * T + 1));
         TE:= 2
     "END" OUT1;
 
     "PROCEDURE" OUT2;
     "IF" T = .6 "THEN"
     OUTPUT(61, "("//, "(" PROBLEM 2")", //,
        "(" NUMBER OF DERIVATIVE CALLS")",
        "("  U(.6, 0)COMPUTED   U(.6, 0)EXACT")", //, 13ZD,
        2(-10Z.7D), "("...")"")", DATA[1] * DATA[8], U[0], EXP(-.09));
 
     I:= 1;
     "FOR" DAT:= 3, 3, 1, 1, "-3, "-6, "-6, 0, 0, 0, 1, .5, 1 / 6 "DO"
     "BEGIN" DATA[I]:= DAT; I:= I + 1 "END";
     T:= 0; Y[1]:= 1; TE:= 1;
     ARK(T, TE, 1, 1, Y, DER1, DATA, OUT1);
     I:= 1;
     "FOR" DAT:= 4, 3, SQRT(8), 500 / 3, DATA[3] / DATA[4], -1, -1,
        0, 0, 0, 1, .5, 1 / 6, 1 / 24 "DO"
     "BEGIN" DATA[I]:= DAT; I:= I + 1 "END";
     M0:= -150; M:= 150; T:= 0; U[0]:= 1;
     "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO"
     U[I]:= U[-I]:= EXP(-(.003 * I) ** 2);
     ARK(T, .6, M0, M, U, DER2, DATA, OUT2)
 "END"
 
     THIS PROGRAM DELIVERS:
 
 
  PROBLEM 1
 
  X  NUMBER OF INTEGRATION STEPS Y(COMPUTED)  Y(EXACT)
 
  1            38                1.7320535    1.7320508...
  2            56                2.2360928    2.2360680...
 
 
  PROBLEM 2
 
  NUMBER OF DERIVATIVE CALLS  U(.6, 0)COMPUTED   U(.6, 0)EXACT
 
            144               .9139326           .9139312...
 
 
1SECTION : 5.2.1.1.1.1.H     (FEBRUARY 1979)                     PAGE 6
 
 
 
 SOURCE TEXT(S):
0"CODE" 33061;
 "PROCEDURE" ARK (T, TE, M0, M, U, DERIVATIVE, DATA, OUT);
 "INTEGER" M0, M;
 "REAL" T, TE;
 "ARRAY" U, DATA;
 "PROCEDURE" DERIVATIVE, OUT;
 
 "BEGIN" "INTEGER" P, N, Q;
     "OWN" "REAL" EC0, EC1, EC2, TAU0, TAU1, TAU2, TAUS, T2;
     "REAL" THETANM1, TAU, BETAN, QINV, ETA;
     "ARRAY" MU, LAMBDA[1:DATA[1]], THETHA[0:DATA[1]], RO, R[M0:M];
     "BOOLEAN" START, STEP1, LAST;
 
     "PROCEDURE" INIVEC(L, U, A, X);               "CODE" 31010;
     "PROCEDURE" MULVEC(L, U, SHIFT, A, B, X);     "CODE" 31020;
     "PROCEDURE" DUPVEC(L, U, SHIFT, A, B);        "CODE" 31030;
     "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
     "PROCEDURE" ELMVEC(L, U, SHIFT, A, B, X);     "CODE" 34020;
     "PROCEDURE" DECSOL(A, N, AUX, B);             "CODE" 34301;
 
     "PROCEDURE" INITIALIZE;
     "BEGIN" "INTEGER" I, J, K, L, N1; "REAL" S, THETA0;
         "ARRAY" ALFA[1:8, 1:DATA[1]+1], TH[1:8], AUX[1:3];
 
         "REAL" "PROCEDURE" LABDA(I, J); "VALUE" I, J; "INTEGER" I, J;
         LABDA:= "IF" P < 3 "THEN" ("IF" J =I-1 "THEN" MUI(I) "ELSE" 0)
                 "ELSE" "IF" P =3 "THEN" ("IF" I =N "THEN" ("IF" J=0
                 "THEN" .25 "ELSE" "IF" J =N - 1 "THEN" .75
                 "ELSE" 0) "ELSE" "IF" J =0 "THEN" ("IF" I =1
                 "THEN" MUI(1) "ELSE" .25) "ELSE" "IF" J =I - 1
                 "THEN" LAMBDA[I] "ELSE" 0) "ELSE" 0;
 
         "REAL" "PROCEDURE" MUI(I); "VALUE" I; "INTEGER" I;
         MUI:= "IF" I =N "THEN" 1 "ELSE"
               "IF" I < 1 ! I > N "THEN" 0 "ELSE"
               "IF" P < 3 "THEN" LAMBDA[I] "ELSE"
               "IF" P =3 "THEN" LAMBDA[I] + .25 "ELSE" 0;
 
         "REAL" "PROCEDURE" SUM(I, A, B, X);
         "VALUE" B; "INTEGER" I, A, B; "REAL" X;
         "BEGIN" "REAL" S; S:= 0;
             "FOR" I:= A "STEP" 1 "UNTIL" B "DO" S:= S + X;
             SUM:= S
         "END" SUM;
                                                               "COMMENT"
1SECTION : 5.2.1.1.1.1.H      (DECEMBER 1979)                     PAGE 7
                                                                  ;
 
 
         N:= DATA[1]; P:= DATA[2]; EC1:= EC2 := 0;
         BETAN:= DATA[3];
         THETANM1:= "IF" P=3 "THEN" .75 "ELSE" 1;
         THETA0:= 1 - THETANM1; S:= 1;
         "FOR" J:= N - 1 "STEP" - 1 "UNTIL" 1 "DO"
         "BEGIN" S:= - S * THETA0 + DATA[N + 10 - J];
             MU[J]:= DATA[N + 11 - J] / S;
             LAMBDA[J]:= MU[J] - THETA0
         "END";
         "FOR" I:= 1 "STEP" 1 "UNTIL" 8 "DO"
         "FOR" J:= 0 "STEP" 1 "UNTIL" N "DO"
         ALFA[I, J + 1]:= "IF" I = 1 "THEN" 1 "ELSE"
           "IF" J = 0 "THEN" 0 "ELSE" "IF" I = 2 ! I = 4 ! I = 8 "THEN"
           MUI(J) ** ENTIER((I + 2) / 3) "ELSE"
           "IF" (I = 3 ! I = 6) & J > 1 "THEN" SUM(L, 1, J-1,
           LABDA(J, L) * MUI(L) ** ENTIER(I / 3)) "ELSE"
           "IF" I = 5 & J > 2 "THEN" SUM(L, 2, J - 1, LABDA(J, L) *
           SUM(K, 1, L - 1, LABDA(L, K) * MUI(K))) "ELSE"
           "IF" I = 7 & J > 1 "THEN" SUM(L, 1, J - 1, LABDA(J, L) *
           MUI(L)) * MUI(J) "ELSE" 0;
         N1:="IF" N < 4 "THEN" N + 1 "ELSE" "IF" N < 7 "THEN" 4
           "ELSE" 8;
         I:= 1;
         "FOR" S:= 1, .5, 1 / 6, 1 / 3, 1 / 24, 1 / 12, .125, .25 "DO"
         "BEGIN" TH[I]:= S; I:= I + 1 "END";
         "IF" P = 3 & N < 7 "THEN" TH[1]:= TH[2]:= 0;
         AUX[2]:= " - 14; DECSOL(ALFA, N1, AUX, TH);
         INIVEC(0, N, THETHA, 0);
         DUPVEC(0, N1 - 1, 1, THETHA, TH);
         "IF" ^ (P = 3 & N < 7) "THEN"
         "BEGIN" THETHA[0]:= THETHA[0] - THETA0;
             THETHA[N - 1]:= THETHA[N - 1] - THETANM1; Q:= P + 1
         "END" "ELSE" Q:= 3;
         QINV:= 1 / Q;
         START:= DATA[8] = 0; DATA[10]:= 0; LAST:= "FALSE";
         DUPVEC(M0, M, 0, R, U); DERIVATIVE(T, R)
     "END" INITIALIZE;
 
1SECTION : 5.2.1.1.1.1.H     (DECEMBER 1975)                     PAGE 8
                                                                  ;
 
 
     "PROCEDURE" LOCAL ERROR CONSTRUCTION(I); "VALUE" I; "INTEGER" I;
     "BEGIN" "IF" THETHA[I] ^= 0 "THEN"
         ELMVEC(M0, M, 0, RO, R, THETHA[I]);
         "IF" I = N "THEN"
         "BEGIN" DATA[9]:= SQRT(VECVEC(M0, M, 0, RO, RO))* TAU;
             EC0:= EC1; EC1:= EC2; EC2:= DATA[9] / TAU ** Q
         "END"
     "END" LEC;
 
     "PROCEDURE" STEPSIZE;
     "BEGIN" "REAL" TAUACC, TAUSTAB, AA, BB, CC, EC;
         ETA:= SQRT(VECVEC(M0, M, 0, U, U)) * DATA[7] + DATA[6];
         "IF" ETA > 0 "THEN"
         "BEGIN" "IF" START "THEN"
             "BEGIN" "IF" DATA[8] = 0 "THEN"
                 "BEGIN" TAUACC:= DATA[5];
                     STEP1:= "TRUE"
                 "END" "ELSE" "IF" STEP1 "THEN"
                 "BEGIN" TAUACC:= (ETA / EC2) ** QINV;
                     "IF" TAUACC > 10 * TAU2 "THEN"
                     TAUACC:= 10 * TAU2 "ELSE" STEP1:= "FALSE"
                 "END" "ELSE"
                 "BEGIN" BB:= (EC2 - EC1) / TAU1; CC:= - BB * T2 + EC2;
                     EC:= BB * T + CC;
                     TAUACC:= "IF" EC < 0 "THEN" TAU2 "ELSE"
                     (ETA / EC) ** QINV;
                     START:= "FALSE"
                 "END"
             "END" "ELSE"
             "BEGIN" AA:= ((EC0 - EC1) / TAU0 + (EC2 - EC1) / TAU1)
                         / (TAU1 + TAU0);
                 BB:= (EC2 - EC1) / TAU1 - (2 * T2 - TAU1) * AA;
                 CC:= - (AA * T2 + BB) * T2 + EC2;
                 EC:= (AA * T + BB) * T + CC;
                 TAUACC:= "IF" EC < 0 "THEN"
                          TAUS "ELSE" (ETA / EC) ** QINV;
                 "IF" TAUACC > 2 * TAUS "THEN" TAUACC:= 2 * TAUS;
                 "IF" TAUACC < TAUS / 2 "THEN" TAUACC:= TAUS / 2
             "END"
         "END" "ELSE" TAUACC:= DATA[5];
          "IF" TAUACC < DATA[5] "THEN" TAUACC:= DATA[5];
         TAUSTAB:= BETAN / DATA[4]; "IF" TAUSTAB < DATA[5] "THEN"
         "BEGIN" DATA[10]:= 1; "GOTO" ENDARK "END";
         TAU:= "IF" TAUACC > TAUSTAB "THEN" TAUSTAB "ELSE" TAUACC;
         TAUS:= TAU; "IF" TAU >= TE - T "THEN"
         "BEGIN" TAU:= TE - T; LAST:= "TRUE" "END";
         TAU0:= TAU1; TAU1:= TAU2; TAU2:= TAU
     "END" STEPSIZE;
 
1SECTION : 5.2.1.1.1.1.H     (DECEMBER 1975)                     PAGE 9
                                                                  ;
 
 
     "PROCEDURE" DIFFERENCE SCHEME;
     "BEGIN" "INTEGER" I, J;
         "REAL" MT, LT;
         MULVEC(M0, M, 0, RO, R, THETHA[0]);
         "IF" P = 3 "THEN" ELMVEC(M0, M, 0, U, R, .25 * TAU);
         "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
         "BEGIN" MT:= MU[I] * TAU; LT:= LAMBDA[I] * TAU;
             "FOR" J:= M0 "STEP" 1 "UNTIL" M "DO"
             R[J]:= LT * R[J] + U[J];
             DERIVATIVE(T + MT, R); LOCAL ERROR CONSTRUCTION(I)
         "END";
         ELMVEC(M0, M, 0, U, R, THETANM1 * TAU);
         DUPVEC(M0, M, 0, R, U); DERIVATIVE(T + TAU, R);
         LOCAL ERROR CONSTRUCTION(N); T2:= T;
         "IF" LAST "THEN"
         "BEGIN" LAST:= "FALSE"; T:= TE "END" "ELSE" T:= T + TAU;
         DATA[8]:= DATA[8]+1
     "END" DIFSCH;
 
     INITIALIZE;
 
   NEXT STEP:
     STEPSIZE; DIFFERENCE SCHEME; OUT;
     "IF" T ^= TE "THEN" "GOTO" NEXT STEP;
 
   ENDARK:
 "END" ARK;
         "EOP"
1SECTION : 5.2.1.1.1.1.I      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE  : EFRK. 
        
        
 AUTHOR: K.DEKKER.
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 740710.
        
        
 BRIEF DESCRIPTION: 
        
     EFRK  SOLVES  AN INITIAL VALUE PROBLEM  FOR A SYSTEM OF FIRST  
     ORDER ORDINARY DIFFERENTIAL EQUATIONS   DU / DT = H(T,U) .     
     EFRK IS A SPECIAL PURPOSE PROCEDURE FOR STIFF EQUATIONS WITH A 
     KNOWN, CLUSTERED EIGENVALUE SPECTRUM.      
        
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEM, 
     SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS, 
     STIFF EQUATION.
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE EFRK READS:   
     "PROCEDURE" EFRK ( T,TE, M0,M, U, SIGMA, PHI, DIAMETER, DERIVATIVE,      
                        K, STEP, R, L, BETA,  THIRDORDER, TOL, OUTPUT );      
     "VALUE" R,L; 
     "INTEGER" M0,M,K,R,L;  
     "REAL" T,TE,SIGMA,PHI,DIAMETER,STEP,TOL;   
     "ARRAY" U,BETA;
     "BOOLEAN" THIRDORDER;  
     "PROCEDURE" DERIVATIVE,OUTPUT;   
     "CODE" 33070;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     T:      <VARIABLE>;    
             THE INDEPENDENT VARIABLE T;
             ENTRY: THE INITIAL VALUE T0;       
             EXIT : THE FINAL VALUE TE; 
     TE:     <ARITHMETIC ETPRESSION>; 
             THE FINAL VALUE OF T (TE>=T);      
     M0:     <ARITHMETIC EXPRESSION>; 
             THE INDEX OF THE FIRST EQUATION;   
     M:      <ARITHMETIC EXPRESSION>; 
             THE INDEX OF THE LAST EQUATION;    
1SECTION : 5.2.1.1.1.1.I      (FEBRUARY 1979)                     PAGE 2      
        
        
        
     U:      <ARRAY IDENTIFIER>;      
             "REAL" "ARRAY" U[M0:M];  
             THE DEPENDENT VARIABLE;  
             ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE  SYSTEM OF      
                    DIFFERENTIAL EQUATIONS AT T = T0;     
             EXIT : THE VALUES OF THE SOLUTION AT T = TE; 
     SIGMA:  <ARITHMETIC EXPRESSION>; 
             THE  MODULUS  OF  THE POINT AT WHICH EXPONENTIAL FITTING IS      
             DESIRED , FOR EXAMPLE AN APPROXIMATION OF THE CENTRE OF THE      
             LEFT HAND CLUSTER;       
     PHI:    <ARITHMETIC EXPRESSION>; 
             THE ARGUMENT OF THE CENTRE OF THE LEFT HAND CLUSTER; IN THE      
             CASE OF  TWO COMPLEX CONJUGATED CLUSTERS , THE  ARGUMENT OF      
             THE CENTRE IN THE SECOND QUADRANT SHOULD BE TAKEN;     
     DIAMETER: <ARITHMETIC EXPRESSION>; 
             THE DIAMETER OF THE LEFT HAND CLUSTER OF EIGENVALUES OF THE      
             JACOBIAN MATRIX  OF THE SYSTEM  OF DIFFERENTIAL  EQUATIONS;      
             IN CASE OF NON-LINEAR EQUATIONS DIAMETER SHOULD HAVE SUCH A      
             VALUE THAT THE VARIATION OF THE EIGENVALUES IN THIS CLUSTER      
             IN THE PERIOD ( T ,T+STEP ) IS LESS THAN HALF THE DIAMETER;      
     DERIVATIVE: <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS: 
             "PROCEDURE" DERIVATIVE(T,U); "REAL" T; "ARRAY" U;      
             THIS PROCEDURE SHOULD DELIVER THE VALUE OF H(T,U) IN THE 
             POINT (T,U) IN THE ARRAY U;
     K:      <VARIABLE>;    
             COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
             FOR EXAMPLE, K MAY BE USED IN THE EXPRESSION FOR TE;   
             ENTRY: AN (ARBIRARY) CHOSEN VALUE K0, E.G. K0=0;       
             EXIT : K0 + THE NUMBER OF INTEGRATION STEPS PERFORMED; 
     STEP:   <ARITHMETIC EXPRESSION>; 
             THE  STEPSIZE  CHOSEN  WILL  BE  AT  MOST  EQUAL TO  STEP ;      
             THIS  STEPSIZE  MAY BE  REDUCED  BY  STABILITY CONSTRAINTS,      
             IMPOSED BY  A  POSITIVE DIAMETER , OR  BY CONSIDERATIONS OF      
             INTERNAL STABILITY (SEE REF[1], PAGE 11);    
     R:      <ARITHMETIC EXPRESSION>; 
             R + L: THE NUMBER OF  EVALUATIONS OF  H(T, U)  ON WHICH THE      
             RUNGE-KUTTA SCHEME IS BASED;       
             FOR R=1,2,>=3 FIRST, SECOND AND THIRD ORDER ACCURACY MAY BE      
             OBTAINED  BY  AN  APPROPRIATE  CHOICE  OF THE  ARRAY  BETA;      
     L:      <ARITHMETIC EXPRESSION>; 
             ENTRY: 
             IF PHI = 4*ARCTAN(1): THE ORDER OF THE EXPONENTIAL FITTING,      
             ELSE TWICE THE ORDER OF THE EXPONENTIAL FITTING;       
             NOTE THAT L SHOULD BE EVEN IN THE LATTER CASE; 
     BETA:   <ARRAY IDENTIFIER>;      
             "REAL" "ARRAY" BETA[0:R+L];
             ENTRY: THE ELEMENTS  BETA[I] , I=0, ... ,R  SHOULD HAVE THE      
                    VALUE OF THE R+1 FIRST COEFFICIENTS OF THE STABILITY      
                    POLYNOMIAL;       
1SECTION : 5.2.1.1.1.1.I      (FEBRUARY 1979)                     PAGE 3      
        
        
        
     THIRDORDER: <BOOLEAN EXPRESSION>;
             IF THIRD ORDER ACCURACY IS DESIRED , THIRDORDER SHOULD HAVE      
             THE VALUE "TRUE" , IN COMBINATION WITH  APPROPRIATE CHOICES      
             OF R (R>=3) AND THE ARRAY BETA ( BETA[I]=1/I!, I=0,1,2,3 );      
             IN ALL OTHER CASES THIRDORDER MUST  HAVE THE VALUE "FALSE";      
     TOL:    <ARITHMETIC EXPRESSION>; 
             AN UPPERBOUND FOR  THE ROUNDING ERRORS  IN THE COMPUTATIONS      
             IN ONE RUNGE-KUTTA STEP ; IN SOME CASES ( E.G. LARGE VALUES      
             OF SIGMA AND R ) TOL WILL CAUSE A DECREASE OF THE STEPSIZE;      
     OUTPUT: <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS; 
             "PROCEDURE" OUTPUT;      
             THIS  PROCEDURE  IS  CALLED  AT THE END OF EACH INTEGRATION      
             STEP ; THE USER CAN ASK FOR OUTPUT OF SOME PARAMETERS , FOR      
             EXAMPLE T, K, U, AND COMPUTE NEW VALUES FOR SIGMA, PHI, AND      
             DIAMETER.      
        
        
 PROCEDURES USED: 
        
     ELMVEC= CP34020,       
     DEC   = CP34300,       
     SOL   = CP34051.       
        
        
 REQUIRED CENTRAL MEMORY : CIRCA 30 + (M - M0) + L * (5 + L).       
        
        
 METHOD AND PERFORMANCE:    
     EFRK IS BASED ON EXPLICIT RUNGE-KUTTA METHODS OF ORDER 1, 2 AND 3,       
     WHICH MAKE USE OF EXPONENTIAL FITTING. AUTOMATIC ERROR CONTROL 
     IS NOT PROVIDED.       
     A DETAILED DESCRIPTION  OF THE METHOD  AND  SOME NUMERICAL EXAMPLES      
     ARE GIVEN IN REF[1]. REF [3], PAGE 170 REPRESENTS A BRIEF SURVEY. A      
     COMPARATIVE TEST  OVER A LARGE CLASS  OF  DIFFERENTIAL EQUATIONS IS      
     GIVEN IN REF [4].      
     FROM THESE RESULTS IT APPEARS  THAT CALLS WITH  THIRDORDER = "TRUE"      
     ARE LESS ADVISABLE.    
        
        
1SECTION : 5.2.1.1.1.1.I      (FEBRUARY 1979)                     PAGE 4      
        
        
        
 REFERENCES:      
        
     [1]. K. DEKKER.
          AN  ALGOL 60  VERSION  OF  EXPONENTIALLY  FITTED   RUNGE-KUTTA      
          METHODS (DUTCH).  
          NR 25 (1972), MATHEMATICAL CENTRE.    
        
     [2]. T. J. DEKKER, P. W. HEMKER AND P. J. VAN DER HOUWEN.      
          COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 1 (DUTCH).
          MC SYLLABUS 15.1, (1972) MATHEMATICAL CENTRE.   
        
     [3]. P. A. BEENTJES, K. DEKKER, H. C. HEMKER, S.P.N. VAN KAMPEN AND      
          G. M. WILLEMS.    
          COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 2 (DUTCH).
          MC SYLLABUS 15.2, (1973) MATHEMATICAL CENTRE.   
        
     [4]. P. A. BEENTJES, K. DEKKER, H. C. HEMKER AND M. V. VELDHUIZEN.       
          COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
          MC SYLLABUS 15.3, (1974) MATHEMATICAL CENTRE.   
        
        
 EXAMPLE OF USE:  
        
     CONSIDER THE SYSTEM OF DIFFERENTIAL EQUATIONS:       
     DY[1]/DX = -Y[1] + Y[1] * Y[2] + .99 * Y[2]
     DY[2]/DX = -1000 * ( -Y[1] + Y[1] * Y[2] + Y[2] )    
     WITH THE INITIAL CONDITIONS AT X = 0:      
     Y[1] = 1 AND Y[2] = 0. (SEE REF[2], PAGE 11).
     THE SOLUTION AT X = 50 IS APPROXIMATELY:   
     Y[1] = .765 878 320 487  AND  Y[2] = .433 710 353  5768.       
     THE FOLLOWING PROGRAM  SHOWS SOME DIFFERENT CALLS  OF THE PROCEDURE      
     EFRK:
        
1SECTION : 5.2.1.1.1.1.I      (AUGUST 1974)                       PAGE 5      
        
        
        
     "BEGIN"      
         "PROCEDURE" EFRK(X,XE,M0,M,Y,SIGMA,PHI,DIAMETER,DERIVATIVE,K,
                                 STEP,R,L,BETA,THIRDORDER,TOL,OUTPUT);
         "CODE" 33070;      
        
         "INTEGER" K,R,L,PASSES;      
         "REAL" X,SIGMA,PHI,TIME,STEP,DIAMETER; 
         "REAL" "ARRAY" Y[1:2],BETA[0:6];       
        
         "PROCEDURE" DER(X,Y); "REAL" X; "ARRAY" Y;       
         "BEGIN" "REAL" Y1,Y2; Y1:=Y[1]; Y2:=Y[2];
             Y[1]:=(Y1+.99)*(Y2-1)+.99; 
             Y[2]:=1000*((1+Y1)*(1-Y2)-1);      
             PASSES:=PASSES+1 
         "END";   
        
         "PROCEDURE" OUT;   
         "BEGIN" "REAL" S;  
             S:=(-1000*Y[1]-1001+Y[2])/2;       
             SIGMA:=ABS(S-SQRT(S*S+10*(Y[2]-1))); 
             DIAMETER:=2*STEP*ABS(1000*(1.99*Y[2]-2*Y[1]*(1-Y[2])));
             "IF" X=50 "THEN" 
             OUTPUT(61,"("4BD,2BD,2(-5ZD),2(4B+.3DB3DB3D),-5ZD.3D,/")",       
             R,L,K,PASSES,Y[1],Y[2],CLOCK-TIME) 
         "END";   
        
         OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
         "("PRINTED BY THIS PROGRAM")",//,      
         "("    THE RESULTS WITH EFRK ARE:")",/,
         "("    R  L     K   DER.EV.        Y[1]             Y[2]")"
         "("         TIME")",/")");   
         PHI:=4*ARCTAN(1); BETA[0]:=BETA[1]:=1; 
         "FOR" R:=1,2,3 "DO" "FOR" L:=1,2,3 "DO"
         "BEGIN" "FOR" K:=2 "STEP" 1 "UNTIL" R "DO"       
             BETA[K]:=BETA[K-1]/K;    
             "FOR" STEP:=1,.1 "DO"    
             "BEGIN" PASSES:=K:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
                 OUT;       
                 EFRK(X,50,1,2,Y,SIGMA,PHI,DIAMETER,DER,K,STEP,R,L,BETA,      
                 R>=3,"-4,OUT);       
             "END"; OUTPUT(61,"("/")"); 
         "END";   
     "END";       
        
1SECTION : 5.2.1.1.1.1.I      (AUGUST 1974)                       PAGE 6      
        
        
        
     THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:   
        
     THE RESULTS WITH EFRK ARE:       
     R  L     K   DER.EV.         Y[1]             Y[2]         TIME
     1  1    237    474     +.765 812 555    +.433 689 306      1.395 
     1  1    501   1002     +.765 847 870    +.433 700 619      3.381 
        
     1  2     52    156     +.765 570 874    +.433 615 119      0.465 
     1  2    501   1503     +.765 848 220    +.433 700 709      4.200 
        
     1  3     52    208     +.765 571 278    +.433 615 202      0.531 
     1  3    500   2000     +.765 848 512    +.433 700 827      4.879 
        
     2  1   3317   9951     +.765 878 320    +.433 710 353     21.808 
     2  1   1050   3150     +.765 878 321    +.433 710 330      7.153 
        
     2  2    174    696     +.765 878 335    +.433 710 335      1.385 
     2  2    501   2004     +.765 878 323    +.433 709 211      4.915 
        
     2  3     57    285     +.765 881 339    +.433 817 185      0.642 
     2  3    501   2505     +.765 878 323    +.433 709 725      5.756 
        
     3  1   7010  28040     +.765 878 320    +.433 710 354     55.298 
     3  1   3255  13020     +.765 878 320    +.433 710 374     25.772 
        
     3  2    949   4745     +.765 878 319    +.433 711 893      8.499 
     3  2   1384   6920     +.765 862 498    +.449 724 830     13.452 
        
     3  3    917   5502     +.765 878 018    +.434 105 184      9.143 
     3  3   1166   6996     +.765 861 696    +.433 705 641     15.512 
        
 SOURCE TEXT(S):  
0"CODE"" 33070;   
 "PROCEDURE" EFRK(T,TE,M0,M,U,SIGMA,PHI,DIAMETER,DERIVATIVE,K,STEP,R,L,       
         BETA,THIRDORDER,TOL,OUTPUT); 
 "VALUE" R,L;     
 "INTEGER" M0,M,K,R,L;      
 "REAL" T,TE,SIGMA,PHI,DIAMETER,STEP,TOL;       
 "ARRAY" U,BETA;  
 "BOOLEAN" THIRDORDER;      
 "PROCEDURE" DERIVATIVE,OUTPUT;       
 "BEGIN" "INTEGER" N;       
     "REAL" THETA0,THETANM1,H,B,B0,PHI0,PHIL,PI,COSPHI,SINPHI,EPS,BETAR;      
     "BOOLEAN" FIRST,LAST,COMPLEX,CHANGE;       
     "INTEGER" "ARRAY" P[1:L];
     "REAL" "ARRAY" MU,LABDA[0:R+L-1],PT[0:R],FAC,BETAC[0:L-1],RL[M0:M],      
                    A[1:L,1:L],AUX[0:3];
     "PROCEDURE" ELMVEC(L,U,SHIFT,A,B,X); "CODE" 34020;   
     "PROCEDURE" SOL(A,N,P,B); "CODE" 34051;    
     "PROCEDURE" DEC(A,N,AUX,P); "CODE" 34300;                 "COMMENT"      
        
1SECTION : 5.2.1.1.1.1.I      (AUGUST 1974)                       PAGE 7      
                                                                  ; 
        
        
     "PROCEDURE" FORM CONSTANTS;      
     "BEGIN" "INTEGER" I;   
         FIRST:="FALSE";    
         FAC[0]:=1; 
         "FOR" I:=1 "STEP" 1 "UNTIL" L-1 "DO" FAC[I]:=I*FAC[I-1];   
         PT[R]:=L*FAC[L-1]; 
         "FOR" I:=1 "STEP" 1 "UNTIL" R "DO"     
         PT[R-I]:=PT[R-I+1]*(L+I)/I   
     "END" FORM CONSTANTS;  
        
     "PROCEDURE" FORM BETA; 
     "BEGIN" "INTEGER" I,J; "REAL" BB,C,D;      
         "IF" FIRST "THEN" FORM CONSTANTS;      
         "IF" L=1 "THEN"    
         "BEGIN" C:=1-EXP(-B);
             "FOR" J:=1 "STEP" 1 "UNTIL" R "DO" C:=BETA[J]-C/B;     
             BETA[R+1]:=C/B 
         "END" "ELSE"       
         "IF" B>40 "THEN"   
         "BEGIN" "FOR" I:=R+1 "STEP" 1 "UNTIL" R+L "DO"   
             "BEGIN" C:=0;  
                 "FOR" J:=0 "STEP" 1 "UNTIL" R "DO"       
                 C:=BETA[J]*PT[J]/(I-J)-C/B;    
                 BETA[I]:=C/B/FAC[L+R-I]/FAC[I-R-1]       
             "END"; 
         "END" "ELSE"       
         "BEGIN" D:=C:=EXP(-B); BETAC[L-1]:=D/FAC[L-1];   
             "FOR" I:=1 "STEP" 1 "UNTIL" L-1 "DO" 
             "BEGIN" C:=B*C/I; D:=D+C; BETAC[L-1-I]:=D/FAC[L-1-I] "END";      
             BB:=1; 
             "FOR" I:=R+1 "STEP" 1 "UNTIL" R+L "DO"       
             "BEGIN" C:=0;  
                 "FOR" J:=0 "STEP" 1 "UNTIL" R "DO"       
                 C:=(BETA[J]-("IF" J<L "THEN" BETAC[J] "ELSE" 0))*  
                    PT[J]/(I-J)-C/B;  
                 BETA[I]:=C/B/FAC[L+R-I]/FAC[I-R-1]+      
                          ("IF" I<L "THEN" BB*BETAC[I] "ELSE" 0);   
                 BB:=BB*B   
             "END"
         "END"    
     "END" FORM BETA;       
        
     "PROCEDURE" SOLUTION OF COMPLEX EQUATIONS; 
     "BEGIN" "INTEGER" I,J,C1,C3;     
         "REAL" C2,E,B1,ZI,COSIPHI,SINIPHI,COSPHIL;       
         "REAL" "ARRAY" D[1:L];                                "COMMENT"      
        
1SECTION : 5.2.1.1.1.1.I      (AUGUST 1974)                       PAGE 8      
                                                                  ; 
        
        
         "PROCEDURE" ELEMENTS OF MATRIX;
         "BEGIN" PHIL:=PHI0;
             COSPHI:=COS(PHIL); SINPHI:=SIN(PHIL);
             COSIPHI:=1; SINIPHI:=0;  
             "FOR" I:=0 "STEP" 1 "UNTIL" L-1 "DO" 
             "BEGIN" C1:=R+1+I; C2:=1;
                 "FOR" J:=L-1 "STEP" -2 "UNTIL" 1 "DO"    
                 "BEGIN" A[J,L-I]:=C2*COSIPHI;  
                     A[J+1,L-I]:=C2*SINIPHI;    
                     C2:=C1*C2; C1:=C1-1
                 "END";     
                 COSPHIL:=COSIPHI*COSPHI-SINIPHI*SINPHI;  
                 SINIPHI:=COSIPHI*SINPHI+SINIPHI*COSPHI;  
                 COSIPHI:=COSPHIL     
             "END"; 
             AUX[2]:=0; DEC(A,L,AUX,P)
         "END" EL OF MAT;   
        
         "PROCEDURE" RIGHTHANDSIDE;   
         "BEGIN" E:=EXP(B*COSPHI);    
             B1:=B*SINPHI-(R+1)*PHIL; 
             COSIPHI:=E*COS(B1); SINIPHI:=E*SIN(B1);      
             B1:=1/B; ZI:=B1**R;      
             "FOR" J:=L "STEP" -2 "UNTIL" 2 "DO"
             "BEGIN" D[J]:=ZI*SINIPHI;
                 D[J-1]:=ZI*COSIPHI;  
                 COSPHIL :=COSIPHI*COSPHI-SINIPHI*SINPHI; 
                 SINIPHI:=COSIPHI*SINPHI+SINIPHI*COSPHI;  
                 COSIPHI:=COSPHIL;    
                 ZI:=ZI*B   
             "END"; 
             COSIPHI:=ZI:=1; SINIPHI:=0;
             "FOR" I:=R "STEP" -1 "UNTIL" 0 "DO"
             "BEGIN" C1:=I; C2:=BETA[I];
                 C3:="IF" 2*I>L-2 "THEN" 2 "ELSE" L-2*I;  
                 COSPHIL :=COSIPHI*COSPHI-SINIPHI*SINPHI; 
                 SINIPHI:=COSIPHI*SINPHI+SINIPHI*COSPHI;  
                 COSIPHI:=COSPHIL;    
                 "FOR" J:=L "STEP" -2 "UNTIL" C3 "DO"     
                 "BEGIN" D[J]:=D[J]+ZI*C2*SINIPHI;
                     D[J-1]:=D[J-1]-ZI*C2*COSIPHI;
                     C2:=C2*C1; C1:=C1-1
                 "END";     
                 ZI:=ZI*B1  
             "END"
         "END" RIGHT HAND SIDE;       
        
         "IF" PHI0^=PHIL "THEN" ELEMENTS OF MATRIX;       
         RIGHTHANDSIDE;     
         SOL(A,L,P,D);      
         "FOR" I:=1 "STEP" 1 "UNTIL" L "DO" BETA[R+I]:=D[L+1-I]*B1  
     "END" SOLOFCOMEQ;                                         "COMMENT"      
        
1SECTION : 5.2.1.1.1.1.I      (AUGUST 1974)                       PAGE 9      
                                                                  ; 
        
        
     "PROCEDURE" COEFFICIENT; 
     "BEGIN" "INTEGER" J,K; "REAL" C; 
         B0:=B; PHI0:=PHI;  
         "IF" B>=1 "THEN"   
         "BEGIN" "IF" COMPLEX "THEN" SOLUTION OF COMPLEX EQUATIONS  
             "ELSE" FORM BETA 
         "END";   
         LABDA[0]:=MU[0]:=0;
         "IF" THIRDORDER "THEN"       
         "BEGIN" THETA0:=.25; THETANM1:=.75;    
             "IF" B<1 "THEN"
             "BEGIN" C:=MU[N-1]:=2/3; LABDA[N-1]:=5/12;   
                 "FOR" J:=N-2 "STEP" -1 "UNTIL" 1 "DO"    
                 "BEGIN" C:=MU[J]:=C/(C-.25)/(N-J+1);     
                     LABDA[J]:=C-.25  
                 "END"      
             "END" "ELSE"   
             "BEGIN" C:=MU[N-1]:=BETA[2]*4/3; LABDA[N-1]:=C-.25;    
                 "FOR" J:=N-2 "STEP" -1 "UNTIL" 1 "DO"    
                 "BEGIN" C:=MU[J]:=C/(C-.25)*BETA[N-J+1]/BETA[N-J]/ 
                            ("IF" J<L "THEN" B "ELSE" 1); 
                     LABDA[J]:=C-.25  
                 "END"      
             "END"
         "END" "ELSE"       
         "BEGIN" THETA0:=0; THETANM1:=1;
             "IF" B<1 "THEN"
             "BEGIN" "FOR" J:=N-1 "STEP" -1 "UNTIL" 1 "DO"
                 MU[J]:=LABDA[J]:=1/(N-J+1)     
             "END" "ELSE"   
             "BEGIN" LABDA[N-1]:=MU[N-1]:=BETA[2];
                 "FOR" J:=N-2 "STEP" -1 "UNTIL" 1 "DO"    
                 MU[J]:=LABDA[J]:=BETA[N-J+1]/BETA[N-J]/  
                        ("IF" J<L "THEN" B "ELSE" 1)      
             "END"
         "END"    
     "END" COEFFICIENT;                                        "COMMENT"      
        
1SECTION : 5.2.1.1.1.1.I      (AUGUST 1974)                      PAGE 10      
                                                                 ;  
        
        
     "PROCEDURE" STEPSIZE;  
     "BEGIN" "REAL" D,HSTAB,HSTABINT; 
         H:=STEP; 
         D:=ABS(SIGMA*SIN(PHI));      
         COMPLEX:=L//2*2=L "AND" 2*D>DIAMETER;  
         "IF" DIAMETER>0 "THEN"       
         HSTAB:=(SIGMA**2/(DIAMETER*(DIAMETER*.25+D)))**(L*.5/R)/   
                BETAR/SIGMA 
         "ELSE" HSTAB:=H;   
         D:= "IF" THIRDORDER "THEN" (2*TOL/EPS/BETA[R])**(1/(N-1))* 
             4**((L-1)/(N-1)) "ELSE" (TOL/EPS)**(1/R)/BETAR;
         HSTABINT:= ABS(D/SIGMA);     
         "IF" H>HSTAB "THEN" H:=HSTAB;
         "IF" H>HSTABINT "THEN" H:=HSTABINT;    
         "IF" T+H>TE*(1-K*EPS) "THEN" 
         "BEGIN" LAST:="TRUE"; H:=TE-T "END";   
         B:=H*SIGMA; D:=DIAMETER*.1*H; D:=D*D;  
         "IF" H<T*EPS  "THEN" "GOTO" ENDOFEFRK; 
         CHANGE:=B0=-1 "OR" ((B-B0)*(B-B0)+B*B0*(PHI-PHI0)*(PHI-PHI0)>D)      
     "END" STEPSIZE;
        
     "PROCEDURE" DIFFERENCESCHEME ;   
     "BEGIN" "INTEGER" I,J; "REAL" MT,LT,THT;   
         I:=-1;   
         NEXTTERM:
         I:=I+1; MT:=MU[I]*H; LT:=LABDA[I]*H;   
         "FOR" J:=M0 "STEP" 1 "UNTIL" M "DO" RL[J]:=U[J]+LT*RL[J];  
         DERIVATIVE(T+MT,RL); 
         "IF" I=0 "OR" I=N-1 "THEN"   
         "BEGIN" THT:="IF" I=0 "THEN" THETA0*H "ELSE" THETANM1*H;   
             ELMVEC(M0,M,0,U,RL,THT)  
         "END";   
         "IF" I<N-1 "THEN" "GOTO" NEXTTERM;     
         T:=T+H   
     "END" DIFFERENCE SCHEME; 
        
     N:=R+L; FIRST:="TRUE"; B0:=-1; BETAR:=BETA[R]**(1/R);
     LAST:="FALSE"; EPS:=2**(-48); PI:=PHI0:=PHIL:=4*ARCTAN(1);     
 NEXTLEVEL:       
     STEPSIZE;    
     "IF" CHANGE "THEN" COEFFICIENT;  
     K:=K+1;      
     DIFFERENCE SCHEME;     
     OUTPUT;      
     "IF" "NOT" LAST "THEN" "GOTO" NEXTLEVEL;   
 ENDOFEFRK:       
 "END" EXPONENTIALLY FITTED RUNGE KUTTA;
         "EOP"    
1SECTION : 5.2.1.1.1.2        (NOVEMBER 1976)                     PAGE 1
  
  
  
 SECTION 5.2.1.1.1.2 CONTAINS SIX ALTERNATIVE PROCEDURES FOR SOLVING
 FIRST-ORDER INITIAL VALUE PROBLEMS WITH THE JACOBIAN MATRIX AVAILABLE.
  
 A.  EFSIRK SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS AN AUTONOMOUS
     SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS DY/DX = F(Y), BY
     MEANS OF AN EXPONENTIALLY FITTED, SEMI-IMPLICIT RUNGE-KUTTA
     METHOD; IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE
     INTEGRATION OF STIFF EQUATIONS.
  
 B.  EFERK SOLVES INITIAL VALUE PROBLEMS, GIVEN AS AN AUTONOMOUS SYSTEM
     OF FIRST ORDER DIFFERENTIAL EQUATIONS, BY MEANS OF AN EXPONENTIALLY
     FITTED, EXPLICIT RUNGE KUTTA METHOD OF THIRD ORDER, WHICH INVOLVES
     THE USE OF THE JACOBIAN MATRIX. AUTOMATIC STEP CONTROL IS PROVIDED.
     IN PARTICULAR  THIS METHOD IS SUITABLE FOR THE INTEGRATION OF STIFF
     DIFFERENTIAL EQUATIONS.
  
 C.  LINIGER1VS SOLVES INITIAL VALUE PROBLEMS, GIVEN AS AN AUTONOMOUS
     SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS, BY MEANS OF AN
     IMPLICIT, FIRST ORDER ACCURATE, EXPONENTIALLY FITTED ONESTEP
     METHOD.
     AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
  
 D.  LINIGER2 SOLVES INITIAL VALUE PROBLEMS, GIVEN AS AN AUTONOMOUS
     SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS, BY MEANS OF AN
     EXPONENTIALLY FITTED ONESTEP METHOD.
     NO AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
  
 E.  GMS SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS AN AUTONOMOUS SYSTEM
     OF FIRST ORDER DIFFERENTIAL EQUATIONS DY / DX = F(Y), BY MEANS OF A
     THIRD ORDER GENERALIZED LINEAR MULTISTEP METHOD.
  
 F.  IMPEX SOLVES AN INITIAL VALUE PROBLEM,GIVEN AS AN AUTONOMOUS SYSTEM
     OF FIRST ORDER DIFFERENTIAL EQUATIONS, BY MEANS OF THE IMPLICIT
     MID-POINT RULE WITH SMOOTHING AND EXTRAPOLATION.
     AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
  
  
 IN PARTICULAR ALL THESE METHODS ARE SUITABLE FOR THE INTEGRATION OF
 STIFF DIFFERENTIAL EQUATIONS.
1SECTION : 5.2.1.1.1.2.A      (AUGUST 1974)                       PAGE 1
  
  
  
 AUTHOR: S.P.N. VAN KAMPEN.
  
  
 INSTITUTE: MATHEMATICAL CENTRE.
  
  
 RECEIVED: 730529.
  
  
 BRIEF DESCRIPTION:
  
     EFSIRK SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS AN AUTONOMOUS
     SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS DY/DX = F(Y), BY
     MEANS OF AN EXPONENTIALLY FITTED, SEMI-IMPLICIT RUNGE-KUTTA
     METHOD; IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE
     INTEGRATION OF STIFF EQUATIONS.
  
  
 KEYWORDS:
  
     DIFFERENTIAL EQUATIONS,
     INITIAL VALUE PROBLEM,
     AUTONOMOUS SYSTEM,
     STIFF EQUATIONS,
     SEMI-IMPLICIT RUNGE-KUTTA METHOD,
     EXPONENTIAL FITTING.
  
  
 CALLING SEQUENCE:
  
     HEADING:
     "PROCEDURE" EFSIRK(X, XE, M, Y, DELTA, DERIVATIVE, JACOBIAN, J,
                        N, AETA, RETA, HMIN, HMAX, LINEAR, OUTPUT);
     "VALUE" M; "INTEGER" M, N;
     "REAL" X, XE, DELTA, AETA, RETA, HMIN, HMAX;
     "PROCEDURE" DERIVATIVE, JACOBIAN, OUTPUT;
     "ARRAY" Y, J;
     "BOOLEAN" LINEAR;
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     X:      <VARIABLE>;
             THE INDEPENDENT VARIABLE X;
             ENTRY: THE INITIAL VALUE X0;
             EXIT : THE END VALUE XE;
     XE:     <ARITHMETIC EXPRESSION>;
             THE END VALUE OF X;
     M:      <ARITHMETIC EXPRESSION>;
             THE NUMBER OF DIFFERENTIAL EQUATIONS;
     Y:      <ARRAY IDENTIFIER>;
             "ARRAY" Y[1 : M];
             THE DEPENDENT VARIABLE;
             DURING THE INTEGRATION PROCESS THE COMPUTED SOLUTION
             AT THE POINT X IS ASSIGNED TO THE ARRAY Y;
             ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM;
1SECTION : 5.2.1.1.1.2.A      (AUGUST 1974)                       PAGE 2
  
  
  
     DELTA:  <ARITHMETIC EXPRESSION>;
             DELTA DENOTES THE REAL PART OF THE POINT AT WHICH
             EXPONENTIAL FITTING IS DESIRED;
             ALTERNATIVES:
             DELTA = (AN ESTIMATE OF) THE REAL PART OF THE, IN  ABSOLUTE
             VALUE, LARGEST EIGENVALUE OF THE JACOBIAN MATRIX OF THE
             SYSTEM;
             DELTA < -10**14, IN ORDER TO OBTAIN ASYMPTOTIC
             STABILITY;
             DELTA = 0, IN ORDER TO OBTAIN A HIGHER ORDER OF ACCURACY IN
             CASE OF LINEAR OR ALMOST LINEAR EQUATIONS;
     DERIVATIVE: <PROCEDURE IDENTIFIER>;
             "PROCEDURE" DERIVATIVE(A); "ARRAY" A;
             WHEN IN EFSIRK DERIVATIVE IS CALLED, A[I] CONTAINS THE
             VALUES OF Y[I];
             UPON COMPLETION OF A CALL OF DERIVATIVE, THE ARRAY A
             SHOULD CONTAIN THE VALUES OF F(Y);
             NOTE THAT THE VARIABLE X SHOULD NOT BE USED IN DERIVATIVE,
             BECAUSE THE DIFFERENTIAL EQUATION IS SUPPOSED TO BE
             AUTONOMOUS;
     JACOBIAN: <PROCEDURE IDENTIFIER>;
             "PROCEDURE" JACOBIAN(J, Y); "ARRAY" J, Y;
             WHEN IN EFSIRK JACOBIAN IS CALLED THE ARRAY Y CONTAINS
             THE VALUES OF THE DEPENDENT VARIABLE;
             UPON COMPLETION OF A CALL OF JACOBIAN THE ARRAY J SHOULD
             CONTAIN THE VALUES OF THE JACOBIAN MATRIX OF F(Y);
     J:      <ARRAY IDENTIFIER>;
             J[1 : M, 1 : M];
             J IS AN AUXILLIARY ARRAY WHICH IS USED IN THE PROCEDURE
             JACOBIAN;
     N:      <VARIABLE>;
             AN INTEGER WHICH COUNTS THE INTEGRATION STEPS;
     AETA, RETA:
             <ARITHMETIC EXPRESSION>;
             REQUIRED ABSOLUTE AND RELATIVE LOCAL ACCURACY;
     HMIN, HMAX:
             <ARITHMETIC EXPRESSION>;
             MINIMAL  AND MAXIMAL  STEPSIZE BY WHICH  THE INTEGRATION IS
             PERFORMED;
     LINEAR: <BOOLEAN EXPRESSION>;
             IF LINEAR = "TRUE" THE PROCEDURE JACOBIAN WILL ONLY BE
             CALLED IF N = 1; THE INTEGRATION WILL THEN BE PERFORMED
             WITH A STEPSIZE HMAX; THE CORRESPONDING REDUCTION
             OF COMPUTING TIME CAN BE EXPLOITED IN CASE OF LINEAR OR
             ALMOST LINEAR EQUATIONS;
     OUTPUT: <PROCEDURE IDENTIFIER>;
             "PROCEDURE" OUTPUT;
             IN OUTPUT ONE MAY PRINT THE VALUES OF E.G. X,
             Y[I], J[K, L] AND N.
  
  
1SECTION : 5.2.1.1.1.2.A      (AUGUST 1974)                       PAGE 3
  
  
  
 DATA AND RESULTS: SEE REF[2] AND [3].
  
  
 PROCEDURES USED:
  
     VECVEC = CP34010,
     MATVEC = CP34011,
     MATMAT = CP34013,
     GSSELM = CP34231,
     SOLELM = CP34061.
  
  
 REQUIRED CENTRAL MEMORY:
  
     EXECUTION FIELD LENGTH: CIRCA M * M + 5 * M.
  
  
 RUNNING TIME:
  
     DEPENDS  STRONGLY  ON  THE  DIFFERENTIAL  EQUATION  TO  BE   SOLVED
  
  
 LANGUAGE: ALGOL 60.
  
  
 METHOD AND PERFORMANCE:
  
     THE PROCEDURE EFSIRK IS AN EXPONENTIALLY FITTED, A-STABLE,
     SEMI-IMPLICIT RUNGE-KUTTA METHOD OF THIRD ORDER (SEE REF[1]
     AND [2]). THE ALGORITHM USES FOR EACH STEP TWO FUNCTION EVALUATIONS
     AND IF LINEAR = "FALSE" ONE EVALUATION OF THE JACOBIAN MATRIX.
     THE STEPSIZE IS NOT DETERMINED BY THE ACCURACY OF THE NUMERICAL
     SOLUTION, BUT BY THE AMOUNT BY WHICH THE GIVEN DIFFERENTIAL
     EQUATION DIFFERS FROM A LINEAR EQUATION (SEE REF[2]).
     THE PROCEDURE DOES NOT REJECT INTEGRATION STEPS.
  
  
 REFERENCES:
  
     [1].P.J. VAN DER HOUWEN.
         ONE-STEP METHODS WITH ADAPTIVE STABILITY FUCNTIONS
         FOR THE INTEGRATION OF DIFFERENTIAL EQUATIONS.
         LECTURE NOTES OF THE CONFERENCE ON
         NUMERISCHE, INSBESONDERE APPROXIMATIONSTHEORETISCHE
         BEHANDLUNG VON FUNKTIONALGLEICHUNGEN.
         OBERWOLFACH, DECEMBER, 3 - 12, 1972.
  
     [2].SYLLABUS COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 2 (DUTCH).
         MATH.CENTR. SYLLABUS 15.2/73.
  
     [3].SYLLABUS COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
         MATH.CENTR. SYLLABUS 15.3/73.
         TO APPEAR IN 1973.
  
  
1SECTION : 5.2.1.1.1.2.A      (AUGUST 1974)                       PAGE 4
  
  
  
 EXAMPLE OF USE:
  
     WE CONSIDER THE DIFFERENTIAL EQUATION
     DY / DX = -EXP(X) * (Y - LN(X)) + 1 / X,
     ON THE INTERVAL  [0.01, 8], WITH  INITIAL VALUE  Y(0.01) = LN(0.01)
     AND ANALYTICAL SOLUTION Y(X) = LN(X);
     FOR THE FIT POINT WE USE THE  EIGENVALUE  OF THE JACOBIAN MATRIX,
     I.E. DELTA = -EXP(X);
  
     "BEGIN"
         "PROCEDURE" EFSIRK(X, XE, M, Y, DELTA, DERIVATIVE, JACOBIAN, J,
                            N, AETA, RETA, HMIN, HMAX, LINEAR, OUTPUT);
         "CODE" 33160;
         "PROCEDURE" DER(Y); "ARRAY" Y;
         "BEGIN" "REAL" Y2; Y2:= Y[2];
             DELTA:= -EXP(Y2); LNX:= LN(Y2);
             Y[1]:= (Y[1] - LNX) * DELTA + 1 / Y2;
             Y[2]:= 1
         "END" DER;
         "PROCEDURE" JAC(J, Y); "ARRAY" J, Y;
         "BEGIN" "REAL" Y2; Y2:= Y[2];
             J[1, 1]:= DELTA;
             J[1, 2]:= (Y[1] - LNX - 1 / Y2) * DELTA - 1 / (Y2 * Y2);
             J[2, 1]:= J[2, 2]:= 0
         "END" JAC;
         "PROCEDURE" OUTP;
         "IF" X = XE "THEN"
         "BEGIN" "REAL" Y1; Y1:= Y[1]; LNX:= LN(X);
             OUTPUT(61, "(""("N = ")", 2ZD,
                        "("   X = ")", +D.D,
                        "("  Y(X) = ")", +D.5D,
                        "("  DELTA = ")", +3ZD.2D, /,
                        "("ABS. ERR. = ")", .2D"+2D,
                        "("  REL. ERR. = ")", .2D"+2D, //")",
                        N, X, Y1, DELTA,
                        ABS(Y1 - LNX), ABS((Y1 - LNX) / LNX));
             "IF" X = 0.4 "THEN" XE:= 8
         "END" OUTP;
         "INTEGER" N;
         "REAL" X, XE, DELTA, LNX;
         "ARRAY" Y[1 : 2], J[1 : 2, 1 : 2];
  
         XE:= 0.4; X:= 0.01; Y[1]:= LN(0.01); Y[2]:= X;
         EFSIRK(X, XE, 2, Y, DELTA, DER, JAC, J,
                N, "-2, "-2, 0.005, 1.5, "FALSE", OUTP)
     "END"
  
     THIS PROGRAM DELIVERS:
  
         N =  10   X = +0.4  Y(X) = -0.91099  DELTA =    -1.44
         ABS. ERR. = .53"-02  REL. ERR. = .58"-02
  
         N = 98    X = +8.0  Y(X) = +2.07911  DELTA = -2980.02
         ABS. ERR. = .33"-03  REL. ERR. = .16"-03
  
  
1SECTION : 5.2.1.1.1.2.A      (AUGUST 1974)                       PAGE 5
  
  
  
 SOURCE TEXT(S):
0"CODE"" 33160;   
     "PROCEDURE" EFSIRK(X, XE, M, Y, DELTA, DERIVATIVE, JACOBIAN, J,
                        N, AETA, RETA, HMIN, HMAX, LINEAR, OUTPUT); 
     "VALUE" M; "INTEGER" M, N;       
     "REAL" X, XE, DELTA, AETA, RETA, HMIN, HMAX;         
     "PROCEDURE" DERIVATIVE, JACOBIAN, OUTPUT;  
     "BOOLEAN" LINEAR;      
     "ARRAY" Y, J;
        
     "BEGIN" "INTEGER" K, L;
         "REAL" STEP, H, MU0, MU1, MU2, THETA0, THETA1, NU1, NU2,   
         NU3, YK, FK, C1, C2, D;      
         "ARRAY" F, K0, LABDA[1 : M], J1[1 : M,1 : M], AUX[1 : 7];  
         "INTEGER" "ARRAY" RI, CI[1 : M];       
         "BOOLEAN" LIN;     
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013; 
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
         "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;       
         "PROCEDURE" SOLELM(A, N, RI, CI, B); "CODE" 34061;         
        
         "REAL" "PROCEDURE" STEPSIZE; 
         "BEGIN" "REAL" DISCR, ETA, S;
             "IF" LINEAR "THEN" S:= H:= HMAX "ELSE"       
             "IF" N = 1 "OR" HMIN = HMAX "THEN" S:= H:= HMIN "ELSE" 
             "BEGIN" ETA:= AETA + RETA * SQRT(VECVEC(1, M, 0, Y, Y));         
                 C1:= NU3 * STEP; "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO"         
                 LABDA[K]:= LABDA[K] + C1 * F[K] - Y[K];  
                 DISCR:= SQRT(VECVEC(1, M, 0, LABDA, LABDA));       
                 S:= H:= (ETA / (0.75 * (ETA + DISCR)) + 0.33) * H; 
                 "IF" H < HMIN "THEN" S:= H:= HMIN "ELSE" 
                 "IF" H > HMAX "THEN" S:= H:= HMAX        
             "END";         
             "IF" X + S > XE "THEN" S:= XE - X; 
             LIN:= STEP = S "AND" LINEAR; STEPSIZE:= S    
         "END" STEPSIZE;    
        
         "PROCEDURE" COEFFICIENT;     
         "BEGIN" "REAL" Z1, E, ALPHA1, A, B;    
             "OWN" "REAL" Z2;         
             Z1:= STEP * DELTA; "IF" N = 1 "THEN" Z2:= Z1 + Z1;     
             "IF" ABS(Z2 - Z1) > " - 6 * ABS(Z1) "OR" Z2 > - 1 "THEN"         
             "BEGIN" A:= Z1 * Z1 + 12; B:= 6 * Z1;        
                 "IF" ABS(Z1) < 0.1 "THEN"      
                 ALPHA1:= (Z1 * Z1 / 140 - 1) * Z1 / 30 "ELSE"      
                 "IF" Z1 < - "14 "THEN" ALPHA1:= 1 / 3 "ELSE"       
                 "IF" Z1 < - 33 "THEN"
                 ALPHA1:= (A + B) / (3 * Z1 * (2 + Z1)) "ELSE"      
                 "BEGIN" E:= "IF" Z1 < 230 "THEN" EXP(Z1) "ELSE" "100;        
                     ALPHA1:= ((A - B) * E - A - B) /     
                              (((2 - Z1) * E - 2 - Z1) * 3 * Z1)    
                 "END";                                        "COMMENT"      
1SECTION : 5.2.1.1.1.2.A      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
                 MU2:= (1 / 3 + ALPHA1) * 0.25; 
                 MU1:= - (1 + ALPHA1) * 0.5;    
                 MU0:= (6 * MU1 + 2) / 9; THETA0:= 0.25;  
                 THETA1:= 0.75; A:= 3 * ALPHA1; 
                 NU3:= (1 + A) / (5 - A) * 0.5; A:= NU3 + NU3;      
                 NU1:= 0.5 - A; NU2:= (1 + A) * 0.75;     
                 Z2:= Z1    
             "END"
         "END" COEFFICIENT; 
        
         "PROCEDURE" DIFFERENCE SCHEME;         
         "BEGIN" DERIVATIVE(F); STEP:= STEPSIZE;
             "IF" "NOT" LINEAR "OR" N = 1 "THEN" JACOBIAN(J, Y);    
             "IF" "NOT" LIN "THEN"    
             "BEGIN" COEFFICIENT;     
                 C1:= STEP * MU1; D:= STEP * STEP * MU2;  
                 "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO"      
                 "BEGIN" "FOR" L:= 1 "STEP" 1 "UNTIL" M "DO"        
                     J1[K,L]:= D * MATMAT(1, M, K, L, J, J) +       
                               C1 * J[K,L];     
                     J1[K,K]:= J1[K,K] + 1      
                 "END";     
                 GSSELM(J1, M, AUX, RI, CI)     
             "END";         
             C1:= STEP * STEP * MU0; D:= STEP * 2 / 3;    
             "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO"
             "BEGIN" K0[K]:= FK:= F[K];         
                 LABDA[K]:= D * FK + C1 * MATVEC(1, M, K, J, F)     
             "END";         
             SOLELM(J1, M, RI, CI, LABDA);      
             "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO" F[K]:= Y[K] + LABDA[K];      
             DERIVATIVE(F); 
             C1:= THETA0 * STEP; C2:= THETA1 * STEP; D:= NU1 * STEP;
             "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO"
             "BEGIN" YK:= Y[K]; FK:= F[K];      
                 LABDA[K]:= YK + D * FK + NU2 * LABDA[K]; 
                 Y[K]:= F[K]:= YK + C1 * K0[K] + C2 * FK  
             "END"
         "END" DIFFERENCE SCHEME;     
        
         AUX[2]:= "-14; AUX[4]:= 8;   
         "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO" F[K]:= Y[K]; 
         N:= 0; OUTPUT; STEP:= 0;     
     NEXT STEP: N:= N + 1;  
         DIFFERENCE SCHEME; X:= X + STEP; OUTPUT;         
         "IF" X < XE "THEN" "GOTO" NEXT STEP    
     "END" EFSIRK;
         "EOP"