1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 1      
        
        
        
 AUTHOR: K.DEKKER.
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 1973/07/31.      
        
        
 BRIEF DESCRIPTION:         
        
     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.
        
        
 KEYWORDS:        
        
     DIFFERENTIAL EQUATIONS,
     INITIAL VALUE PROBLEMS,
     STIFF EQUATIONS,       
     EXPONENTIAL FITTING,   
     EXPLICIT RUNGE KUTTA METHODS.    
        
        
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE EFERK READS:  
     "PROCEDURE" EFERK(X,XE,M,Y,SIGMA,PHI,DERIVATIVE,J,JACOBIAN,    
                   K,L,AUT,AETA,RETA,HMIN,HMAX,LINEAR,OUTPUT);      
     "VALUE" L;   
     "INTEGER" M,K,L;       
     "REAL" X,XE,SIGMA,PHI,AETA,RETA,HMIN,HMAX; 
     "ARRAY" Y,J; 
     "BOOLEAN" AUT,LINEAR;  
     "PROCEDURE" DERIVATIVE,JACOBIAN,OUTPUT;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:      <VARIABLE>;    
             THE INDEPENDENT VARIABLE X;        
             CAN BE USED IN DERIVATIVE, JACOBIAN, OUTPUT, ETC.;     
             ENTRY: THE INITIAL VALUE X0;       
             EXIT : THE FINAL VALUE XE;         
     XE:     <ARITHMETIC EXPRESSION>; 
             THE FINAL VALUE OF X (XE>=X);      
     M:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF EQUATIONS; 
     Y:      <ARRAY IDENTIFIER>;      
             "REAL" "ARRAY" Y[1:M];   
             THE DEPENDENT VARIABLE;  
             ENTRY: THE  INITIAL  VALUES  OF  THE SYSTEM OF DIFFERENTIAL      
                    EQUATIONS: Y[I] AT X=X0;    
             EXIT : THE FINAL VALUES OF THE SOLUTION: Y[I] AT X=XE; 
     SIGMA:  <ARITHMETIC EXPRESSION>; 
             THE  MODULUS  OF  THE POINT AT WHICH EXPONENTIAL FITTING IS      
             DESIRED, FOR EXAMPLE THE LARGEST NEGATIVE EIGENVALUE OF THE      
             JACOBIAN MATRIX OF THE SYSTEM OF DIFFERENTIAL EQUATIONS;         
     PHI:    <ARITHMETIC EXPRESSION>; 
             THE  ARGUMENT  OF THE  COMPLEX POINT  AT WHICH  EXPONENTIAL      
             FITTING IS DESIRED;      
     DERIVATIVE: <PROCEDURE IDENTIFIER>;        
             THE HEADING OF THIS PROCEDURE READS:         
             "PROCEDURE" DERIVATIVE(Y); "ARRAY" Y;        
             THIS  PROCEDURE  SHOULD  DELIVER THE RIGHT HAND SIDE OF THE      
             I-TH DIFFERENTIAL EQUATION AT THE POINT (Y) AS Y[I];   
     J:      <ARRAY IDENTIFIER>;      
             "REAL" "ARRAY" J[1:M,1:M];         
             THE JACOBIAN MATRIX OF THE SYSTEM; 
             THE ARRAY J SHOULD BE UPDATED IN THE PROCEDURE JACOBIAN;         
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 3      
        
        
        
     JACOBIAN: <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:         
             "PROCEDURE" JACOBIAN(J,Y); "ARRAY" J,Y;      
             IN THIS PROCEDURE  THE JACOBIAN AT THE POINT (Y)  HAS TO BE      
             ASSIGNED TO THE ARRAY J; 
     K:      <VARIABLE>;    
             COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
             FOR EXAMPLE, MAY BE USED IN THE EXPRESSION FOR XE;     
     L:      <ARITHMETIC EXPRESSION>; 
             ENTRY:         
             IF PHI = 4*ARCTAN(1): THE ORDER OF THE EXPONENTIAL FITTING,      
             ELSE TWICE THE ORDER OF THE EXPONENTIAL FITTING;       
     AUT:    <BOOLEAN EXPRESSION>;    
             IF THE SYSTEM HAS BEEN WRITTEN IN AUTONOMOUS FORM BY ADDING      
             THE EQUATION DY[M]/DX = 1 TO THE SYSTEM , THEN AUT MAY HAVE      
             THE VALUE "FALSE", ELSE AUT SHOULD HAVE THE VALUE "TRUE";        
     AETA:   <ARITHMETIC EXPRESSION>; 
             REQUIRED  ABSOLUTE  PRECISION  IN THE  INTEGRATION PROCESS;      
             AETA HAS TO BE POSITIVE; 
     RETA:   <ARITHMETIC EXPRESSION>; 
             REQUIRED  RELATIVE  PRECISION  IN THE  INTEGRATION PROCESS;      
             RETA HAS TO BE POSITIVE; 
     HMIN:   <ARITHMETIC EXPRESSION>; 
             THE STEPLENGTH CHOSEN WILL BE AT LEAST EQUAL TO HMIN;  
     HMAX:   <ARITHMETIC EXPRESSION>; 
             THE STEPLENGTH CHOSEN WILL BE AT MOST EQUAL TO HMAX;   
     LINEAR: <ARITHMETIC EXPRESSION>; 
             THE PROCEDURE JACOBIAN IS CALLED  ONLY IF LINEAR="FALSE" OR      
             K=0; SO IF THE SYSTEM IS LINEAR , LINEAR MAY HAVE THE VALUE      
             "TRUE";        
     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, K, Y.         
        
        
 DATA AND RESULTS: SEE EXAMPLE OF USE, AND REF[4].        
        
        
 PROCEDURES USED: 
        
     VECVEC= CP34010,       
     MATVEC= CP34011,       
     DEC   = CP34300,       
     SOL   = CP34051.       
        
        
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 4      
        
        
        
 REQUIRED CENTRAL MEMORY:   
     EXECUTION FIELD LENGTH: CIRCA 30 + 4 * M + L * (5+L).
        
        
 RUNNING TIME: DEPENDS  STRONGLY  ON THE DIFFERENTIAL EQUATION TO SOLVE.      
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE:    
        
    THE PROCEDURE EFERK IS AN EXPONENTIALLY FITTED, SEMI-EXPLICIT RUNGE       
    KUTTA METHOD OF THIRD ORDER ( SEE REF [1] AND [3] ) . THE ALGORITHM       
    USES FOR EACH STEP TWO FUNCTION EVALUATIONS AND IF LINEAR = "FALSE"       
    ONE EVALUATION OF THE JACOBIAN MATRIX.THE STEPSIZE IS DETERMINED BY       
    AN ESTIMATION  OF THE LOCAL TRUNCATION ERROR  BASED ON THE RESIDUAL       
    FUNCTION (SEE REF[3]). INTEGRATION STEPS ARE NOT REJECTED.      
        
        
 REFERENCES:      
        
    [1]. P.J.VAN DER HOUWEN.
         ONE-STEP  METHODS  WITH ADAPTIVE  STABILITY FUNCTIONS  FOR THE       
         INTEGRATION OF DIFFERENTIAL EQUATIONS. 
         LECTURES NOTES OF THE CONFERENCE ON NUMERISCHE, INSBESONDERE         
         APPROXIMATIONSTHEORETISCHE BEHANDLUNG VON FUNKTIONAL-      
         GLEICHUNGEN.       
         OBERWOLFACH, DECEMBER, 3 -12, 1972.    
        
    [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]. (TO APPEAR).       
          COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).        
          MC SYLLABUS 15.3, (1973) 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      
     EFERK,AND ILLUSTRATES THE ACCURACIES WHICH MAY BE OBTAINED BY THEM:      
        
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 5      
        
        
        
     "BEGIN"      
         "PROCEDURE" EFERK(X,XE,M,Y,SIGMA,PHI,DERIVATIVE,J,JACOBIAN,
                       K,L,AUT,AETA,RETA,HMIN,HMAX,LINEAR,OUTPUT);  
         "CODE" 33120;      
        
         "INTEGER" K,PASSES,PASJAC;   
         "REAL" X,SIGMA,PHI,TIME,TOL; 
         "REAL" "ARRAY" Y[1:2],J[1:2,1:2];      
        
         "PROCEDURE" DER(Y); "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" JACOBIAN(J,Y); "ARRAY" J,Y;
         "BEGIN" J[1,1]:=Y[2]-1; J[1,2]:=.99+Y[1];        
             J[2,1]:=1000*(1-Y[2]); J[2,2]:=-1000*(1+Y[1]);         
             SIGMA:=ABS(J[2,2]+J[1,1]-SQRT((J[2,2]-J[1,1])**2+      
                    4*J[2,1]*J[1,2]))/2;        
             PASJAC:=PASJAC+1         
         "END" JACOBIAN;    
        
         "PROCEDURE" OUT;   
         "IF" X=50 "THEN"   
        OUTPUT(61,"("3(-5ZD),2(4B+.3DB3DB3D),-5ZD.3D,/")",K,PASSES, 
         PASJAC,Y[1],Y[2],CLOCK-TIME);
        
         OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
         "("PRINTED BY THIS PROGRAM")",//,      
         "("    THE RESULTS WITH EFERK -FIRST ORDER FIT- ARE:")",/, 
         "("     K   DER.EV. JAC.EV.       Y[1]             Y[2]")" 
         "("         TIME")",/")");   
         PHI:=4*ARCTAN(1);  
         "FOR" TOL:=1,"-1,"-2,"-3 "DO"
         "BEGIN" PASSES:=PASJAC:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;         
             EFERK(X,50,2,Y,SIGMA,PHI,DER,J,JACOBIAN,K,1,"TRUE",TOL,
                                            TOL,"-6,50,"FALSE",OUT);
         "END";   
     "END";       
        
     THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:   
        
     THE RESULTS WITH EFERK -FIRST ORDER FIT- ARE:        
      K   DER.EV. JAC.EV.       Y[1]             Y[2]         TIME  
     93    186     93    +.765 883 211    +.428 752 781      1.170  
    105    210    105    +.765 878 445    +.433 569 561      1.296  
    147    294    147    +.765 878 317    +.433 708 489      1.834  
    266    532    266    +.765 878 320    +.433 710 229      3.297  
        
        
        
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 6      
        
        
        
 SOURCE TEXT(S):  
0"CODE"" 33120;   
 "PROCEDURE" EFERK(X,XE,M,Y,SIGMA,PHI,DERIVATIVE,J,JACOBIAN,        
               K,L,AUT,AETA,RETA,HMIN,HMAX,LINEAR,OUTPUT);
 "VALUE" L; "INTEGER" M,K,L;
 "REAL" X,XE,SIGMA,PHI,AETA,RETA,HMIN,HMAX; "ARRAY" Y,J;  
 "BOOLEAN" AUT,LINEAR; "PROCEDURE" DERIVATIVE,JACOBIAN,OUTPUT;      
 "BEGIN" "INTEGER" M1,I;    
     "REAL" H,B,B0,PHI0,COSPHI,SINPHI,ETA,DISCR,FAC,PI;   
     "BOOLEAN" CHANGE,LAST; 
     "INTEGER" "ARRAY" P[1:L];        
     "REAL" "ARRAY" BETA,BETHA[0:L],BETAC[0:L+3],K0,D,D1,D2[1:M],   
         A[1:L,1:L],AUX[1:3];         
     "REAL" "PROCEDURE" VECVEC(L,U,SHIFT,A,B); "CODE" 34010;        
     "REAL" "PROCEDURE" MATVEC(L,U,I,A,B); "CODE" 34011;  
     "PROCEDURE" DEC(A,N,AUX,P); "CODE" 34300;  
     "PROCEDURE" SOL(A,N,P,B); "CODE" 34051;    
     "REAL" "PROCEDURE" SUM(I,L,U,T); "VALUE" L,U; "INTEGER" I,L,U; 
     "REAL" T;    
     "BEGIN" "REAL" S; S:=0;
         "FOR" I:=L "STEP" 1 "UNTIL" U "DO" S:=S+T;       
         SUM:=S   
     "END";       
     "PROCEDURE" FORMBETA;  
     "IF" L=1 "THEN"        
     "BEGIN" BETHA[1]:=(.5-(1-(1-EXP(-B))/B)/B)/B;        
         BETA[1]:=(1/6-BETHA[1])/B    
     "END" "ELSE" 
     "IF" L=2 "THEN"        
     "BEGIN" "REAL" E,EMIN1; E:=EXP(-B); EMIN1:=E-1;      
         BETHA[1]:=(1-(3+E+4*EMIN1/B)/B)/B;     
         BETHA[2]:=(.5-(2+E+3*EMIN1/B)/B)/B/B;  
         BETA[2]:=(1/6-BETHA[1])/B/B; 
         BETA[1]:=(1/3-(1.5-(4+E+5*EMIN1/B)/B)/B)/B       
     "END" "ELSE" 
     "BEGIN" "REAL" B0,B1,B2,A0,A1,A2,A3,C,D;   
         BETAC[L-1]:=C:=D:=EXP(-B)/FAC;         
         "FOR" I:=L-1 "STEP" -1 "UNTIL" 1 "DO"  
         "BEGIN" C:=I*B*C/(L-I); BETAC[I-1]:=D:=D*I+C "END";        
         B2:=.5-BETAC[2];   
         B1:=(1-BETAC[1])*(L+1)/B;    
         B0:=(1-BETAC[0])*(L+2)*(L+1)*.5/B/B;   
         A3:=1/6-BETAC[3];  
         A2:=B2*(L+1)/B;    
         A1:=B1*(L+2)*.5/B; 
         A0:=B0*(L+3)/3/B;  
         D:=L/B;  
         "FOR" I:=1 "STEP" 1 "UNTIL" L "DO"     
         "BEGIN"BETA[I]:=(A3/I-A2/(I+1)+A1/(I+2)-A0/(I+3))*D+BETAC[I+3];      
             BETHA[I]:=(B2/I-B1/(I+1)+B0/(I+2))*D+BETAC[I+2];       
             D:=D*(L-I)/I/B;
         "END"    
     "END" FORMBETA;                                           "COMMENT"      
        
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 7      
                                                                  ; 
        
        
     "PROCEDURE" SOLUTIONOFCOMPLEXEQUATIONS;    
     "IF" L=2 "THEN"        
     "BEGIN" "REAL" COS2PHI,COSA,SINA,E,ZI;     
         PHI0:=PHI; COSPHI:=COS(PHI0); SINPHI:=SIN(PHI0); 
         E:=EXP(B*COSPHI); ZI:=B*SINPHI-3*PHI0; 
         SINA:=("IF" ABS(SINPHI)<"-6 "THEN" -E*(B+3)      
               "ELSE" E*SIN(ZI)/SINPHI);        
         COS2PHI:=2*COSPHI*COSPHI-1;  
         BETHA[2]:=(.5+(2*COSPHI+(1+2*COS2PHI+SINA)/B)/B)/B/B;      
         SINA:=("IF" ABS(SINPHI)<"-6 "THEN" E*(B+4)       
               "ELSE" SINA*COSPHI-E*COS(ZI));   
         BETHA[1]:=-(COSPHI+(1+2*COS2PHI+(4*COSPHI*COS2PHI+SINA)    
                   /B)/B)/B;
         BETA[1]:=BETHA[2]+2*COSPHI*(BETHA[1]-1/6)/B;     
         BETA[2]:=(1/6-BETHA[1])/B/B  
     "END" "ELSE" 
        
     "BEGIN" "INTEGER" J,C1;
         "REAL" C2,E,ZI,COSIPHI,SINIPHI,COSPHIL;
         "REAL" "ARRAY" D[1:L];       
         "PROCEDURE" ELEMENTS OF MATRIX;        
         "BEGIN" PHI0:=PHI; 
             COSPHI:=COS(PHI0); SINPHI:=SIN(PHI0);        
             COSIPHI:=1; SINIPHI:=0;  
             "FOR" I:=0 "STEP" 1 "UNTIL" L-1 "DO"         
             "BEGIN" C1:=4+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:=C2*C1; 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" RIGHT HAND SIDE; 
         "BEGIN" E:=EXP(B*COSPHI);    
             ZI:=B*SINPHI-4*PHI0;     
             COSIPHI:=E*COS(ZI); SINIPHI:=E*SIN(ZI);      
             ZI:=1/B/B/B;   
             "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";                                            "COMMENT"      
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 8      
                                                                  ; 
        
        
             SINIPHI:=2*SINPHI*COSPHI;
             COSIPHI:=2*COSPHI*COSPHI-1;        
             COSPHIL:=COSPHI*(2*COSIPHI-1);     
             D[L]:=D[L]+SINPHI*(1/6+(COSPHI+(1+2*COSIPHI*(1+2*COSPHI/B))      
                         /B)/B);      
             D[L-1]:=D[L-1]-COSPHI/6-(.5*COSIPHI+(COSPHIL+(2*COSIPHI*         
                         COSIPHI-1)/B)/B)/B;    
             D[L-2]:=D[L-2]+SINPHI*(.5+(2*COSPHI+(2*COSIPHI+1)/B)/B);         
             D[L-3]:=D[L-3]-.5*COSPHI-(COSIPHI+COSPHIL/B)/B;        
             "IF" L<5 "THEN" "GOTO" END;        
             D[L-4]:=D[L-4]+SINPHI+SINIPHI/B;   
             D[L-5]:=D[L-5]-COSPHI-COSIPHI/B;   
             "IF" L<7 "THEN" "GOTO" END;        
             D[L-6]:=D[L-6]+SINPHI;   
             D[L-7]:=D[L-7]-COSPHI;   
         END:     
         "END" RHS;         
         "IF" PHI0^=PHI "THEN" ELEMENTS OF MATRIX;        
         RIGHT HAND SIDE;   
         SOL(A,L,P,D);      
         ZI:=1/B; 
         "FOR" I:=1 "STEP" 1 "UNTIL" L "DO"     
         "BEGIN" BETA[I]:=D[L+1-I]*ZI;
             BETHA[I]:=(I+3)*BETA[I]; 
             ZI:=ZI/B       
         "END"    
     "END" SOLOFEQCOM;      
        
     "PROCEDURE" COEFFICIENT;         
     "BEGIN" B0:=B:=ABS(H*SIGMA);     
         "IF" B>=.1 "THEN"  
         "BEGIN" "IF" PHI^=PI "AND" L=2 "OR" ABS(PHI-PI)>.01 "THEN" 
             SOLUTION OF COMPLEX EQUATIONS "ELSE" FORMBETA
         "END" "ELSE"       
         "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" L "DO"       
             "BEGIN" BETHA[I]:=BETA[I-1];       
                 BETA[I]:=BETA[I-1]/(I+3);      
             "END"
         "END"    
     "END" COEFFICIENT;     
        
     "PROCEDURE" LOCAL ERROR BOUND;   
     ETA:=AETA+RETA*SQRT(VECVEC(1,M1,0,Y,Y));   
        
     "PROCEDURE" STEPSIZE;  
     "BEGIN" LOCAL ERROR BOUND;       
         "IF" K=0 "THEN"    
         "BEGIN" DISCR:=SQRT(VECVEC(1,M1,0,D,D)); H:=ETA/DISCR      
         "END" "ELSE"       
         "BEGIN" DISCR:=H*SQRT(SUM(I,1,M1,(D[I]-D2[I])**2))/ETA;    
             H:=H*("IF" LINEAR "THEN" 4/(4+DISCR)+.5      
                               "ELSE" 4/(3+DISCR)+1/3)    
         "END";                                                "COMMENT"      
1SECTION : 5.2.1.1.1.2.B      (AUGUST 1974)                       PAGE 9      
                                                                  ; 
        
        
         "IF" H<HMIN "THEN" H:=HMIN;  
         "IF" H>HMAX "THEN" H:=HMAX;  
         B:=ABS(H*SIGMA);   
         CHANGE:=ABS(1-B/B0)>.05 "OR" PHI^=PHI0;
         "IF" 1.1*H>=XE-X "THEN"      
         "BEGIN" CHANGE:=LAST:="TRUE"; H:=XE-X "END";     
         "IF" "NOT" CHANGE "THEN" H:=H*B0/B     
     "END" STEPSIZE;        
        
     "PROCEDURE" DIFFERENCE SCHEME;   
     "BEGIN" "INTEGER" K;   
         "REAL" BETAI,BETHAI;         
         "IF" M1<M "THEN"   
         "BEGIN" D2[M]:=1; K0[M]:=Y[M]+2*H/3; Y[M]:=Y[M]+.25*H "END";         
         "FOR" K:=1 "STEP" 1 "UNTIL" M1 "DO"    
         "BEGIN" K0[K]:=Y[K]+2*H/3*D[K];        
             Y[K]:=Y[K]+.25*H*D[K];   
             D1[K]:=H*MATVEC(1,M,K,J,D);        
             D2[K]:=D1[K]+D[K]        
         "END";   
         "FOR" I:=0 "STEP" 1 "UNTIL" L "DO"     
         "BEGIN" BETAI:=4*BETA[I]/3; BETHAI:=BETHA[I];    
             "FOR" K:=1 "STEP" 1 "UNTIL" M1 "DO" D[K]:=H*D1[K];     
             "FOR" K:=1 "STEP" 1 "UNTIL" M1 "DO"
             "BEGIN" K0[K]:=K0[K]+BETAI*D[K];   
                 D1[K]:=MATVEC(1,M1,K,J,D);     
                 D2[K]:=D2[K]+BETHAI*D1[K]      
             "END"
         "END";   
         DERIVATIVE(K0);    
         "FOR" K:=1 "STEP" 1 "UNTIL" M "DO" Y[K]:=Y[K]+.75*H*K0[K]  
     "END" DIFF SCHEME;     
        
     B0:=PHI0:=-1; PI:=4*ARCTAN(1);   
     BETAC[L]:=BETAC[L+1]:=BETAC[L+2]:=BETAC[L+3]:=0;     
     BETA[0]:=1/6; BETHA[0]:=.5;      
     FAC:=1; "FOR" I:=2 "STEP" 1 "UNTIL" L-1 "DO" FAC:=I*FAC;       
     M1:= "IF" AUT "THEN" M "ELSE" M-1;         
     K:=0; LAST:="FALSE";   
  NEXT LEVEL:     
     "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" D[I]:=Y[I];       
     DERIVATIVE(D);         
     "IF" "NOT" LINEAR "OR" K=0 "THEN" JACOBIAN(J,Y);     
     STEPSIZE;    
     "IF" CHANGE "THEN" COEFFICIENT;  
     OUTPUT;      
     DIFFERENCE SCHEME;     
     K:=K+1;      
     X:=X+H;      
     "IF" "NOT" LAST "THEN" "GOTO" NEXT LEVEL;  
 END OF EFERK: OUTPUT;      
 "END" EFERK;     
         "EOP"    
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 1      
        
        
        
 AUTHOR: K.DEKKER.
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 1973/07/16.      
        
        
 BRIEF DESCRIPTION:         
        
     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. 
     IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION  OF STIFF      
     DIFFERENTIAL EQUATIONS.
        
        
 KEYWORDS:        
        
     DIFFERENTIAL EQUATIONS,
     INITIAL VALUE PROBLEMS,
     STIFF EQUATIONS,       
     EXPONENTIAL FITTING,   
     IMPLICIT ONESTEP METHODS.        
        
        
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 2      
        
        
        
 CALLING SEQUENCES:         
        
     THE HEADING OF THE PROCEDURE LINIGER2 READS:         
     "PROCEDURE" LINIGER2(X,XE,M,Y,SIGMA1,SIGMA2,F,EVALUATE,J,      
                    JACOBIAN,K,ITMAX,STEP,AETA,RETA,OUTPUT);        
     "INTEGER" M,K,ITMAX;   
     "REAL" X,XE,SIGMA1,SIGMA2,STEP,AETA,RETA;  
     "ARRAY" Y,J; 
     "BOOLEAN" "PROCEDURE" EVALUATE;  
     "REAL" "PROCEDURE" F;  
     "PROCEDURE" JACOBIAN,OUTPUT;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:      <VARIABLE>;    
             THE INDEPENDENT VARIABLE X;        
             CAN BE USED IN F, JACOBIAN, OUTPUT, ETC.;    
             ENTRY: THE INITIAL VALUE X0;       
             EXIT : THE FINAL VALUE XE;         
     XE:     <ARITHMETIC EXPRESSION>; 
             THE FINAL VALUE OF X (XE>=X);      
     M:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF EQUATIONS; 
     Y:      <ARRAY IDENTIFIER>;      
             "ARRAY" Y[1:M];
             THE DEPENDENT VARIABLE;  
             ENTRY: THE  INITIAL  VALUES  OF  THE SYSTEM OF DIFFERENTIAL      
                    EQUATIONS: Y[I] AT X=X0;    
             EXIT : THE FINAL VALUES OF THE SOLUTION: Y[I] AT X=XE; 
     SIGMA1: <ARITHMETIC EXPRESSION>; 
             THE  MODULUS  OF  THE POINT AT WHICH EXPONENTIAL FITTING IS      
             DESIRED;  THIS  POINT  MAY BE COMPLEX OR REAL AND NEGATIVE;      
     SIGMA2: <ARITHMETIC EXPRESSION>; 
             SIGMA2  MAY DEFINE  THREE  DIFFERENT TYPES  OF  EXPONENTIAL      
             FITTING: FITTING IN TWO COMPLEX CONJUGATED POINTS , FITTING      
             IN  TWO  REAL  NEGATIVE  POINTS , OR  FITTING  IN ONE POINT      
             COMBINED WITH THIRD ORDER ACCURACY;
             IF THIRD ORDER ACCURACY IS DESIRED , SIGMA2 SHOULD HAVE THE      
             VALUE 0;       
             IF FITTING  IN A SECOND NEGATIVE POINT  IS DESIRED , SIGMA2      
             SHOULD HAVE THE VALUE OF THE MODULUS OF THIS POINT;    
             IF FITTING IN  TWO COMPLEX CONJUGATED  POINTS  IS DESIRED ,      
             THEN  SIGMA2 SHOULD BE  MINUS THE VALUE  OF THE ARGUMENT OF      
             THE POINT IN THE SECOND QUADRANT (THUS A VALUE BETWEEN - PI      
             AND - PI/2);   
     F:      <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS:         
             "REAL" "PROCEDURE" F(I); "INTEGER" I;        
             THIS  PROCEDURE  SHOULD  DELIVER THE RIGHT HAND SIDE OF THE      
             I-TH DIFFERENTIAL EQUATION AS F;   
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 3      
        
        
        
     EVALUATE: <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:         
             "BOOLEAN" "PROCEDURE" EVALUATE(ITNUM); "INTEGER" ITNUM;
             EVALUATE  SHOULD  HAVE  THE  VALUE "TRUE", IF IT IS DESIRED      
             THAT  THE JACOBIAN OF THE SYSTEM IS UPDATED IN THE ITNUM-TH      
             ITERATION STEP OF THE NEWTON PROCESS;        
      J:     <ARRAY IDENTIFIER>;      
             "ARRAY" J[1:M,1:M];      
             THE JACOBIAN MATRIX OF THE SYSTEM; 
             THE ARRAY J SHOULD BE UPDATED IN THE PROCEDURE JACOBIAN;         
     JACOBIAN: <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:         
             "PROCEDURE" JACOBIAN(J,Y); "ARRAY" J,Y;      
             IN THIS PROCEDURE  THE JACOBIAN  HAS TO BE ASSIGNED  TO THE      
             THE ARRAY J , OR AN APPROXIMATION OF THE JACOBIAN , IF ONLY      
             SECOND ORDER ACCURACY IS REQUIRED; 
     K:      <VARIABLE>;    
             COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
             FOR EXAMPLE, CAN BE USED IN EVALUATE;        
     ITMAX:  <ARITHMETIC EXPRESSION>; 
             AN  UPPERBOUND  FOR  THE NUMBER  OF ITERATIONS  IN NEWTON'S      
             PROCESS, USED TO SOLVE THE IMPLICIT EQUATIONS;         
     STEP:   <ARITHMETIC EXPRESSION>; 
             THE LENGTH OF THE INTEGRATION STEP, TO BE PRESCRIBED BY THE      
             THE USER;      
     AETA:   <ARITHMETIC EXPRESSION>; 
             REQUIRED ABSOLUTE PRECISION  IN THE NEWTON PROCESS, USED TO      
             SOLVE THE IMPLICIT EQUATIONS;      
     RETA:   <ARITHMETIC EXPRESSION>; 
             REQUIRED RELATIVE PRECISION  IN THE NEWTON PROCESS, USED TO      
             SOLVE THE IMPLICIT EQUATIONS;      
     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 THE PARAMETERS , FOR      
             EXAMPLE X, K, Y.         
        
        
 DATA AND RESULTS: SEE EXAMPLE OF USE, AND REF[3].        
        
        
 PROCEDURES USED: 
        
     VECVEC= CP34010,       
     MATVEC= CP34011,       
     MATMAT= CP34013,       
     DEC   = CP34300,       
     SOL   = CP34051.       
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: CIRCA 20 + M * (4+M).        
        
        
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 4      
        
        
        
 RUNNING TIME: DEPENDS  STRONGLY  ON THE DIFFERENTIAL EQUATION TO SOLVE.      
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE:    
        
     LINIGER2: INTEGRATES  THE  SYSTEM OF DIFFERENTIAL EQUATIONS FROM X0      
               UNTIL XE, BY MEANS OF A SECOND ORDER FORMULA (IF SIGMA2=0      
               AND EVALUATE="TRUE" EVEN THIRD ORDER).     
     SEE ALSO REF[1] AND REF[2].      
        
 REFERENCES:      
        
     [1]. W.LINIGER AND R.A.WILLOUGHBY.         
          EFFICIENT  INTEGRATION  METHODS  FOR STIFF SYSTEMS OF ORDINARY      
          DIFFERENTIAL EQUATIONS.     
          SIAM J. NUM. ANAL. 7 (1970) PAGE 47.  
        
     [2]. T.J.DEKKER, P.W.HEMKER AND P.W.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.   
        
        
 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 2487  AND  Y[2] = .433 710 353  5768.      
     THE FOLLOWING PROGRAM  SHOWS SOME DIFFERENT CALLS OF  THE PROCEDURE      
     LINIGER2, AND ILLUSTRATES THE ACCURACY WHICH MAY BE OBTAINED BY IT:      
        
     "BEGIN"      
         "PROCEDURE" LINIGER2(X,XE,M,Y,SIGMA1,SIGMA2,F,EVALUATE,J,  
                        JACOBIAN,K,ITMAX,STEP,AETA,RETA,OUTPUT);    
         "CODE" 33131;      
        
         "INTEGER" K,ITMAX,PASSES,PASJAC;       
         "REAL" X,SIGMA,STEP,TIME;    
         "REAL" "ARRAY" Y[1:2],J[1:2,1:2];      
        
         "REAL" "PROCEDURE" F(I); "INTEGER" I;  
         "IF" I=1 "THEN" F:=(Y[1]+.99)*(Y[2]-1)+.99 "ELSE"
         "BEGIN" PASSES:=PASSES+1; F:=1000*((1+Y[1])*(1-Y[2])-1) "END";       
        
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 5      
        
        
        
         "PROCEDURE" JACOBIAN(J,Y); "ARRAY" J,Y;
         "BEGIN" J[1,1]:=Y[2]-1; J[1,2]:=.99+Y[1];        
             J[2,1]:=1000*(1-Y[2]); J[2,2]:=-1000*(1+Y[1]);         
             SIGMA:=ABS(J[2,2]+J[1,1]-SQRT((J[2,2]-J[1,1])**2+      
                    4*J[2,1]*J[1,2]))/2;        
             PASJAC:=PASJAC+1         
         "END" JACOBIAN;    
         "BOOLEAN" "PROCEDURE" EVALUATE1(I); "INTEGER" I; 
         EVALUATE1:= I=1;   
         "BOOLEAN" "PROCEDURE" EVALUATE2(I); "INTEGER" I; 
         EVALUATE2:= "TRUE";
         "PROCEDURE" OUT;   
         "IF" X=50 "THEN"   
         OUTPUT(61,"("3(-4ZDB),2(4B+.3DB3DB3D),-5ZD.3D,/")",K,PASSES,         
         PASJAC,Y[1],Y[2],CLOCK-TIME);
         OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
         "("PRINTED BY THIS PROGRAM")",//,      
         "("    THE RESULTS WITH LINIGER2 -SECOND ORDER- ARE:")",/, 
         "("     K   DER.EV. JAC.EV.       Y[1]             Y[2]")" 
         "("         TIME")",/")");   
         "FOR" STEP:=10,1 "DO"        
         "FOR" ITMAX:=1,3 "DO"        
         "BEGIN" PASSES:=PASJAC:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;         
             LINIGER2(X,50,2,Y,SIGMA,0,F,EVALUATE1,J,JACOBIAN,K,ITMAX,        
                                                    STEP,"-4,"-4,OUT);        
         "END";   
         OUTPUT(61,"("//,   
         "("    THE RESULTS WITH LINIGER2 -THIRD ORDER- ARE:")",/,  
         "("     K   DER.EV. JAC.EV.       Y[1]             Y[2]")" 
         "("         TIME")",/")");   
         "FOR" STEP:=10,1 "DO"        
         "FOR" ITMAX:=1,3 "DO"        
         "BEGIN" PASSES:=PASJAC:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;         
             LINIGER2(X,50,2,Y,SIGMA,0,F,EVALUATE2,J,JACOBIAN,K,ITMAX,        
                                                    STEP,"-4,"-4,OUT);        
         "END";   
     "END";       
        
     THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:   
        
     THE RESULTS WITH LINIGER2 -SECOND ORDER- ARE:        
      K   DER.EV. JAC.EV.       Y[1]             Y[2]         TIME  
      5      5      5     +.766 392 210    +.434 218 863      0.092 
      5     15      5     +.765 755 853    +.433 671 223      0.175 
     50     50     50     +.765 884 310    +.433 715 687      0.949 
     50    101     50     +.765 877 388    +.433 710 059      1.494 
        
        
     THE RESULTS WITH LINIGER2 -THIRD ORDER- ARE:         
      K   DER.EV. JAC.EV.       Y[1]             Y[2]         TIME  
      5      5      5     +.766 392 210    +.434 218 863      0.092 
      5     15     15     +.765 882 250    +.433 711 614      0.300 
     50     50     50     +.765 884 310    +.433 715 687      0.949 
     50    101    101     +.765 878 873    +.433 710 531      2.080 
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 6      
        
        
        
 SOURCE TEXT(S):  
        
0"CODE"" 33131;   
 "PROCEDURE" LINIGER2(X,XE,M,Y,SIGMA1,SIGMA2,F,EVALUATE,J,
                JACOBIAN,K,ITMAX,STEP,AETA,RETA,OUTPUT);  
 "INTEGER" M,K,ITMAX;       
 "REAL" X,XE,SIGMA1,SIGMA2,STEP,AETA,RETA;      
 "ARRAY" Y,J;     
 "BOOLEAN" "PROCEDURE" EVALUATE;      
 "REAL" "PROCEDURE" F;      
 "PROCEDURE" JACOBIAN,OUTPUT;         
        
 "BEGIN" "INTEGER" I;       
     "REAL" H,HL,B1,B2,P,Q,C0,C1,C2,C3,C4;      
     "BOOLEAN" LAST;        
     "INTEGER" "ARRAY" PI[1:M];       
     "REAL" "ARRAY" DY,YL,FL[1:M],A[1:M,1:M],AUX[1:3];    
     "REAL" "PROCEDURE" VECVEC(L,U,SHIFT,A,B); "CODE" 34010;        
     "REAL" "PROCEDURE" MATVEC(L,U,I,A,B); "CODE" 34011;  
     "REAL" "PROCEDURE" MATMAT(L,U,I,J,A,B); "CODE" 34013;
     "PROCEDURE" DEC(A,N,AUX,P); "CODE" 34300;  
     "PROCEDURE" SOL(A,N,P,B); "CODE" 34051;    
        
     "PROCEDURE" STEPSIZE;  
     "BEGIN" H:=STEP;       
         "IF" 1.1*H>=XE-X "THEN"      
         "BEGIN" LAST:="TRUE"; H:=XE-X; X:=XE   
         "END" "ELSE" X:=X+H
     "END" STEPSIZE;        
        
     "PROCEDURE" COEFFICIENT;         
     "BEGIN" "REAL" R1,R2,EX,ZETA,ETA,SINL,COSL,SINH,COSH,D;        
         "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;    
         "IF" X>40 "THEN" R:=X/(X-2) "ELSE"     
         "BEGIN" EX:=EXP(-X); R:=X*(1-EX)/(X-2+(X+2)*EX) "END";     
        
         B1:=H*SIGMA1;      
         B2:=H*SIGMA2;      
         "IF" B1<.1 "THEN" "BEGIN" P:=0; Q:=1/3; "GOTO" OUT "END";  
         "IF" B2<0 "THEN" "GOTO" COMPLEX;       
         "IF" B1<1 "OR" B2<.1 "THEN" "GOTO" THIRDORDER;   
         "IF" ABS(B1-B2)<B1*B1*"-6 "THEN" "GOTO" DOUBLEFIT;         
        
         R1:=R(B1)*B1; R2:=R(B2)*B2;  
         D:=B2*R1-B1*R2;    
         P:=2*(R2-R1)/D;    
         Q:=2*(B2-B1)/D;    
         "GOTO" OUT;                                           "COMMENT"      
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 7      
                                                                  ; 
        
        
 THIRDORDER: Q:=1/3;        
         P:=R(B1)/3-2/B1;   
         "GOTO" OUT;        
 DOUBLEFIT: B1:=.5*(B1+B2); 
         R1:=R(B1);         
         "IF" B1>40 "THEN" EX:=0;     
         R2:=B1/(1-EX); R2:=1-EX*R2*R2;         
         Q:=1/(R1*R1*R2);   
         P:=R1*Q-2/B1;      
         "GOTO" OUT;        
 COMPLEX: ETA:=ABS(B1*SIN(SIGMA2));   
         ZETA:=ABS(B1*COS(SIGMA2));   
         "IF" ETA<B1*B1*"-6 "THEN"    
         "BEGIN" B1:=B2:=ZETA; "GOTO" DOUBLEFIT "END";    
         "IF" ZETA>40 "THEN"
         "BEGIN" P:=1-4*ZETA/B1/B1; Q:=4*(1-ZETA)/B1/B1+1 "END" "ELSE"        
         "BEGIN" EX:=EXP(ZETA);       
             SINL:=SIN(ETA); COSL:=COS(ETA);    
             SINH:=.5*(EX-1/EX); COSH:=.5*(EX+1/EX);      
             D:=ETA*(COSH-COSL)-.5*B1*B1*SINL;  
             P:=(ZETA*SINL+ETA*SINH-4*ZETA*ETA/B1/B1*(COSH-COSL))/D;
             Q:=ETA*((COSH-COSL-ZETA*SINH-ETA*SINL)*4/B1/B1+COSH+COSL)/D      
         "END";   
 OUT:    C0:=.25*H*H*(P+Q); 
         C1:=.5*H*(1+P);    
         C2:=H-C1;
         C3:=.25*H*H*(Q-P); 
         C4:=.5*H*P;        
         ELEMENTS OF MATRIX 
     "END" COEFFICIENT;     
        
     "PROCEDURE" ELEMENTS OF MATRIX;  
     "BEGIN" "INTEGER" K;   
         "FOR" I:=1 "STEP" 1 "UNTIL" M "DO"     
         "BEGIN" "FOR" K:=1 "STEP" 1 "UNTIL" M "DO"       
             A[I,K]:=C0*MATMAT(1,M,I,K,J,J)-C1*J[I,K];    
             A[I,I]:=A[I,I]+1         
         "END";   
         AUX[2]:=0; DEC(A,M,AUX,PI)   
     "END" ELOFMAT;                                            "COMMENT"      
1SECTION : 5.2.1.1.1.2.D      (AUGUST 1974)                       PAGE 8      
                                                                  ; 
        
        
     "PROCEDURE" NEWTON ITERATION;    
     "BEGIN" "INTEGER" ITNUM; "REAL" JFL,ETA,DISCR;       
         ITNUM:=0;
 NEXT:   ITNUM:=ITNUM+1;    
         "IF" EVALUATE(ITNUM) "THEN"  
         "BEGIN" JACOBIAN(J,Y); COEFFICIENT "END"         
         "ELSE" "IF" ITNUM=1 "AND" H^=HL "THEN" COEFFICIENT;        
         "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" FL[I]:=F(I);  
         "IF" ITNUM=1 "THEN"
         "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" M "DO"       
             "BEGIN" JFL:=MATVEC(1,M,I,J,FL);   
                 DY[I]:=H*(FL[I]-C4*JFL);       
                 YL[I]:=Y[I]+C2*FL[I]+C3*JFL    
             "END"
         "END" "ELSE"       
         "FOR" I:=1 "STEP" 1 "UNTIL" M "DO"     
         DY[I]:=YL[I]-Y[I]+C1*FL[I]-C0*MATVEC(1,M,I,J,FL);
         SOL(A,M,PI,DY);    
         "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" Y[I]:=Y[I]+DY[I];       
         "IF" ITNUM<ITMAX "THEN"      
         "BEGIN" ETA:=SQRT(VECVEC(1,M,0,Y,Y))*RETA+AETA;  
             DISCR:=SQRT(VECVEC(1,M,0,DY,DY));  
             "IF" ETA<DISCR "THEN" "GOTO" NEXT  
         "END"    
     "END" NEWTON;
        
     LAST:="FALSE"; K:=0; HL:=0;      
 NEXT LEVEL:      
     K:=K+1;      
     STEPSIZE;    
     NEWTON ITERATION;      
     HL:=H;       
     OUTPUT;      
     "IF" "NOT" LAST "THEN" "GOTO" NEXT LEVEL   
 "END" LINIGER2;  
         "EOP"    
1SECTION : 5.2.1.1.1.3        (AUGUST 1974)                       PAGE 1      
        
        
        
 SECTION 5.2.1.1.1.3 CONTAINS TWO  ALTERNATIVE PROCEDURES FOR SOLVING         
 FIRST-ORDER INITIAL VALUE PROBLEMS WITH SEVERAL DERIVATIVES AVAILABLE.       
        
 A.  MODIFIED TAYLOR SOLVES AN INITIAL ( BOUNDARY ) VALUE PROBLEM, GIVEN      
     AS A SYSTEM OF  FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF A       
     ONE-STEP TAYLOR-METHOD.
     IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF LARGE       
     SYSTEMS ARISING FROM PARTIAL DIFFERENTIAL EQUATIONS, PROVIDED THAT       
     HIGHER ORDER DERIVATIVES CAN BE EASILY OBTAINED.     
        
 B.  EXPONENTIALLY FITTED TAYLOR  SOLVES AN INITIAL VALUE PROBLEM, GIVEN      
     AS A SYSTEM  OF  FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF A      
     ONE-STEP  TAYLOR-METHOD . AUTOMATIC  STEPSIZE  CONTROL IS PROVIDED.      
     IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF  STIFF      
     DIFFERENTIAL EQUATIONS , PROVIDED THAT HIGHER ORDER DERIVATIVES CAN      
     BE EASILY OBTAINED.    
1SECTION : 5.2.1.1.1.3.A      (AUGUST 1974)                       PAGE 1      
        
        
        
 AUTHORS: P.J. VAN DER HOUWEN AND P.A.BEENTJES. 
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730616.
        
        
 BRIEF DESCRIPTION:         
        
     MODIFIED TAYLOR SOLVES AN INITIAL ( BOUNDARY ) VALUE PROBLEM, GIVEN      
     AS A SYSTEM  OF  FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF A      
     ONE-STEP TAYLOR-METHOD.
     IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF  LARGE      
     SYSTEMS ARISING FROM PARTIAL DIFFERENTIAL EQUATIONS , PROVIDED THAT      
     HIGHER ORDER DERIVATIVES CAN BE EASILY OBTAINED.     
        
        
 KEYWORDS:        
        
     DIFFERENTIAL EQUATIONS,
     INITIAL (BOUNDARY) VALUE PROBLEMS,         
     ONE-STEP TAYLOR-METHOD.
        
        
1SECTION : 5.2.1.1.1.3.A      (DECEMBER 1975)                     PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE" MODIFIED TAYLOR (T,TE,M0,M,U,SIGMA,TAUMIN,I,DERIVATIVE,      
                                K,DATA,ALFA,NORM,AETA,RETA,ETA,RHO,OUT);      
     "INTEGER" M0,M,I,K,NORM;         
     "REAL" T,TE,SIGMA,TAUMIN,ALFA,AETA,RETA,RHO;         
     "ARRAY" U,DATA;        
     "PROCEDURE" DERIVATIVE,OUT;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     T:      <VARIABLE>;    
             THE INDEPENDENT VARIABLE T;        
             MAY BE USED IN DERIVATIVE, SIGMA ETC.;       
             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 TO  BE      
             SOLVED;        
     U:      <ARRAY IDENTIFIER>;      
             "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 SPECTRAL RADIUS OF THE JACOBIAN MATRIX WITH RESPECT
             TO THOSE EIGENVALUES WHICH ARE LOCATED IN THE LEFT     
             HALFPLANE;     
             IF  SIGMA  TENDS  TO  INFINITY , PROCEDURE  MODIFIED TAYLOR      
             TERMINATES;    
     TAUMIN: <ARITHMETIC EXPRESSION>; 
             MINIMAL STEP LENGTH BY WHICH THE INTEGRATION IS  PERFORMED;      
             HOWEVER,ACTUAL STEPSIZES WILL ALWAYS BE WITHIN THE INTERVAL      
             [MIN(HMIN,HSTAB),HSTAB],WHERE HSTAB(= DATA[0]/SIGMA) IS THE      
             STEPLENGTH PRESCRIBED BY STABILITY CONSIDERATIONS;     
     I:      <VARIABLE>;    
             A JENSEN PARAMETER FOR PROCEDURE DERIVATIVE; 
             MAY BE USED IN M0 AND M; 
     DERIVATIVE: <PROCEDURE IDENTIFIER>;        
             THE HEADING OF THIS PROCEDURE READS:         
             "PROCEDURE" DERIVATIVE(I,A); "INTEGER" I; "ARRAY" A;   
             WHEN THIS PROCEDURE IS CALLED, ARRAY A[M0 : M] CONTAINS THE      
             COMPONENTS OF THE (I-1)-ST DERIVATIVE OF U  AT THE POINT T;      
             UPON COMPLETION OF DERIVATIVE, ARRAY A SHOULD  CONTAIN  THE      
             COMPONENTS  OF THE I-TH  DERIVATIVE  OF U AT  THE POINT  T;      
     K:      <VARIABLE>;    
             INDICATES  THE  NUMBER  OF  INTEGRATION  STEPS   PERFORMED;      
             ENTRY: K = 0;  
1SECTION : 5.2.1.1.1.3.A      (AUGUST 1974)                       PAGE 3      
        
        
        
     DATA:   <ARRAY IDENTIFIER>;      
             "ARRAY" DATA[-2 : DATA[-2]];       
             ENTRY:         
             DATA[-2]: THE  ORDER OF THE  HIGHEST DERIVATIVE  UPON WHICH      
                       THE TAYLOR METHOD IS BASED;        
             DATA[-1]: ORDER OF ACCURACY OF THE METHOD;   
             DATA[0] : STABILITY PARAMETER;     
             DATA[1] , ... , DATA[DATA[-2]] :  POLYNOMIAL  COEFFICIENTS;      
             FOR FURTHER EXPLANATION AND POSSIBLE VALUES OF THE ELEMENTS      
             OF ARRAY DATA SEE REFERENCES [2] AND [3];    
     ALFA:   <ARITHMETIC EXPRESSION>; 
             GROWTH FACTOR FOR THE INTEGRATION STEP LENGTH;         
     NORM:   <ARITHMETIC EXPRESSION>; 
             IF NORM = 1 DISCREPANCY AND TOLERANCE ARE ESTIMATED IN  THE      
             MAXIMUM NORM, OTHERWISE IN THE EUCLIDIAN NORM;         
     AETA,RETA: <ARITHMETIC EXPRESSION>;        
             DESIRED ABSOLUTE AND RELATIVE ACCURACY;      
             IF BOTH  AETA AND RETA  ARE NEGATIVE , ACCURACY  CONDITIONS      
             WILL BE IGNORED;         
     ETA,RHO: <VARIABLE>;   
             COMPUTED TOLERANCE AND DISCREPANCY;
     OUT:    <PROCEDURE IDENTIFIER>;  
             THE  HEADING  OF  THIS  PROCEDURE  READS : "PROCEDURE" OUT;      
             THROUGH  THIS  PROCEDURE  THE VALUES AFTER EACH INTEGRATION      
             STEP  OF  FOR INSTANCE  T, U, ETA AND RHO  ARE  ACCESSIBLE.      
        
        
 DATA AND RESULTS:
        
     FOR FURTHER EXPLANATION OF THE PARAMETERS AETA, RETA, ETA, RHO, M0,      
     M AND THE ARRAY DATA SEE REFERENCES [2] AND [3].     
     AS  FOR  THE  INDICES  M0 AND M THE FOLLOWING MAY BE REMARKED: 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 MODIFIED TAYLOR , THIS  MAY  BE  REALIZED  BY  INTEGER      
     PROCEDURES M0 AND M WHICH  ARE  DEFINED  AS  FUNCTIONS OF I, K  AND      
     DATA[-2].    
        
        
 PROCEDURES USED: VECVEC = CP34010.   
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: CIRCA 75 + M - M0. 
        
        
 RUNNING TIME:    
        
     DEPENDS  STRONGLY  ON  THE  DIFFERENTIAL  EQUATION  TO  BE  SOLVED.      
        
        
 LANGUAGE: ALGOL 60.        
        
        
1SECTION : 5.2.1.1.1.3.A      (AUGUST 1974)                       PAGE 4      
        
        
        
 METHOD AND PERFORMANCE: SEE REFERENCES.        
        
 REFERENCES:      
        
     [1].P.J. VAN DER HOUWEN.         
         ONE-STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS I,      
         POLYNOMIAL METHODS,TW REPORT 119,      
         MATHEMATICAL CENTRE, AMSTERDAM (1970). 
        
     [2].P.J. VAN DER HOUWEN, P.BEENTJES, K.DEKKER AND E.SLAGT      
         ONE-STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS III,    
         NUMERICAL EXAMPLES, TW REPORT 130/71,  
         MATHEMATICAL CENTRE, AMSTERDAM (1971). 
        
     [3].P.J. VAN DER HOUWEN, J.KOK.  
         NUMERICAL SOLUTION OF A MINIMAX PROBLEM, TW REPORT 123/71, 
         MATHEMATICAL CENTRE, AMSTERDAM (1971). 
        
        
 EXAMPLE OF USE:  
        
     THE  SOLUTION AT T=EXP(1) AND T=EXP(2) OF THE DIFFERENTIAL EQUATION      
     DU/DT=-EXP(T)*(U-LN(T)) + 1/T WITH INITIAL CONDITION U(.01)=LN(.01)      
     AND  ANALYTICAL SOLUTION U(T) = LN(T), MAY BE OBTAINED  AS FOLLOWS:      
        
     "BEGIN" "INTEGER" I,K;"REAL" T,TE,ETA,RHO,EXPT,LNT,C0,C1,C2,C3;
         "ARRAY" U[0:0],DATA[-2:4];   
         "PROCEDURE" OP;"IF" T=TE "THEN"        
         OUTPUT(61,"(""("NUMBER OF STEPS:")",3ZD,/,       
                      "("SOLUTION:  T= ")",+D.5D,         
                      "("    U(T) = ")",+D.7D,//")",K,T,U[0]);      
         "PROCEDURE" DER(I,A);"INTEGER" I;"ARRAY" A;      
         "BEGIN" "IF" I=1 "THEN"      
             "BEGIN" EXPT:=EXP(T);LNT:=LN(T);C0:=A[0];    
                 C1:=A[0]:=-EXPT*C0+1/T+EXPT*LNT
             "END";         
             "IF" I=2 "THEN" C2:=A[0]:=EXPT*(LNT+1/T-C0-C1)-1/T/T;  
             "IF" I=3 "THEN" C3:=A[0]:=         
             EXPT*(LNT+2/T-C0-2*C1-C2-1/T/T)+2/T/T/T;     
             "IF" I=4 "THEN" A[0]:=C3-2*(1+3/T)/T/T/T+    
             EXPT*((1-(2-2/T)/T)/T-C1-C2*2-C3)  
         "END";   
         "PROCEDURE" MODIFIED TAYLOR(T,TE,M0,M,U,SIGMA,TAUMIN,I,    
                     DERIVATIVE,K,DATA,ALFA,NORM,AETA,RETA,ETA,RHO,OUT);      
         "CODE" 33040;      
         I:=-2;"FOR" T:=4,3,6.025,1,.5,1/6,.018455702 "DO"
         "BEGIN" DATA[I]:=T;I:=I+1 "END";       
         T:=U[0]:="-2;K:=0;"FOR" TE:=EXP(1),TE*TE "DO"    
         MODIFIED TAYLOR(T,TE,0,0,U,EXP(T),"-4,I,DER,K,DATA,1.5,1,"-5,        
                        "-4,ETA,RHO,OP)         
     "END"        
        
1SECTION : 5.2.1.1.1.3.A      (AUGUST 1974)                       PAGE 5      
        
        
        
     THIS PROGRAM DELIVERS: 
        
     NUMBER OF STEPS:  46   
     SOLUTION:  T= +2.71828    U(T) = +1.0000285
        
     NUMBER OF STEPS: 424   
     SOLUTION:  T= +7.38906    U(T) = +1.9999967
        
        
 SOURCE TEXT(S):  
0"CODE"" 33040;   
 "PROCEDURE" MODIFIED TAYLOR(T,TE,M0,M,U,SIGMA,TAUMIN,I,DERIVATIVE,K,         
                   DATA,ALFA,NORM,AETA,RETA,ETA,RHO,OUT); 
 "INTEGER" M0,M,I,K,NORM;   
 "REAL" T,TE,SIGMA,TAUMIN,ALFA,AETA,RETA,ETA,RHO;         
 "ARRAY" U,DATA;  
 "PROCEDURE" DERIVATIVE,OUT;
        
 "BEGIN" I:=0;    
     "BEGIN" "INTEGER" N,P,Q;         
         "OWN" "REAL" EC0,EC1,EC2,TAU0,TAU1,TAU2,TAUS,T2; 
         "REAL" T0,TAU,TAUI,TAUEC,ECL,BETAN,GAMMA;        
         "REAL" "ARRAY" C[M0:M],BETA,BETHA[1:DATA[-2]];   
         "BOOLEAN" START,STEP1,LAST;  
         "REAL" "PROCEDURE" VECVEC(L,U,SHIFT,A,B); "CODE" 34010;    
        
         "PROCEDURE" COEFFICIENT;     
         "BEGIN" "INTEGER" J;"REAL" IFAC;       
             IFAC:=1; GAMMA:=.5; N:=DATA[-2]; P:=DATA[-1];
             BETAN:=DATA[0]; Q:= "IF" P<N "THEN" P+1 "ELSE" N;      
             "FOR" J:=1 "STEP" 1 "UNTIL" N "DO" 
             "BEGIN" BETA[J]:=DATA[J]; IFAC:=IFAC/J;      
                 BETHA[J]:=IFAC-BETA[J]         
             "END";         
             "IF" P=N "THEN" BETHA[N]:=IFAC     
         "END";                                                "COMMENT"      
        
1SECTION : 5.2.1.1.1.3.A      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
        
         "REAL" "PROCEDURE" NORMFUNCTION(NORM,W);         
         "INTEGER" NORM; "ARRAY" W;   
         "BEGIN" "INTEGER" J; "REAL" S,X;       
             S:=0;
             "IF" NORM=1 "THEN"       
             "BEGIN" "FOR" J:=M0 "STEP" 1 "UNTIL" M "DO"  
                 "BEGIN"  X:=ABS(W[J]); "IF" X>S "THEN" S:=X "END"  
             "END" "ELSE"   
             S:=SQRT(VECVEC(M0,M,0,W,W));       
             NORMFUNCTION:=S
         "END";   
        
         "PROCEDURE" LOCAL ERROR BOUND;         
         ETA:=AETA+RETA * NORMFUNCTION(NORM,U); 
        
         "PROCEDURE" LOCAL ERROR CONSTRUCTION(I);"INTEGER" I;       
         "BEGIN" "IF" I=P "THEN" "BEGIN" ECL:=0;TAUEC:=1 "END";     
             "IF" I>P+1 "THEN" TAUEC:=TAUEC*TAU;
             ECL:=ECL+ABS(BETHA[I])*TAUEC*NORMFUNCTION(NORM,C);     
             "IF" I=N "THEN"
             "BEGIN" EC0:=EC1;EC1:=EC2;EC2:=ECL;
                 RHO:=ECL*TAU**Q      
             "END"
         "END";   
        
         "PROCEDURE" STEPSIZE;        
         "BEGIN" "REAL" TAUACC,TAUSTAB,AA,BB,CC,EC;       
             LOCAL ERROR BOUND;       
             "IF" ETA>0 "THEN"        
             "BEGIN" "IF" START "THEN"
                 "BEGIN" "IF" K=0 "THEN"        
                     "BEGIN" "INTEGER" J;       
                         "FOR" J:=M0 "STEP" 1 "UNTIL" M "DO" C[J]:=U[J];      
                         I:=1; DERIVATIVE(I,C); 
                         TAUACC:=ETA/NORMFUNCTION(NORM,C);
                         STEP1:="TRUE"
                     "END" "ELSE"     
                     "IF" STEP1 "THEN"
                     "BEGIN" TAUACC:=(ETA/RHO)**(1/Q)*TAU2;         
                         "IF" TAUACC>10*TAU2 "THEN"       
                         TAUACC:=10*TAU2 "ELSE" STEP1:="FALSE"      
                     "END" "ELSE"     
                     "BEGIN" BB:=(EC2-EC1)/TAU1; CC:=EC2-BB*T2;     
                         EC:=BB*T+CC; 
                         TAUACC:="IF" EC<0 "THEN" TAU2 "ELSE"       
                         (ETA/EC)**(1/Q);       
                         START:="FALSE"         
                     "END"  
1SECTION : 5.2.1.1.1.3.A      (AUGUST 1974)                       PAGE 7      
        
        
        
                 "END" "ELSE"         
                 "BEGIN" AA:=((EC0-EC1)/TAU0+(EC2-EC1)/TAU1)/       
                             (TAU1+TAU0);       
                     BB:=(EC2-EC1)/TAU1-AA*(2*T2-TAU1);   
                     CC:=EC2-T2*(BB+AA*T2); EC:=CC+T*(BB+T*AA);     
                     TAUACC:="IF" EC<0 "THEN" TAUS        
                             "ELSE" (ETA/EC)**(1/Q);      
                     "IF" TAUACC>ALFA*TAUS "THEN" TAUACC:=ALFA*TAUS;
                     "IF" TAUACC<GAMMA*TAUS "THEN" TAUACC:=GAMMA*TAUS;        
                 "END"      
             "END" "ELSE" TAUACC:=TE-T;         
             "IF" TAUACC<TAUMIN "THEN" TAUACC:=TAUMIN;    
             TAUSTAB:=BETAN/SIGMA;    
             "IF" TAUSTAB<"-12 * (T-T0) "THEN"  
             "BEGIN" OUT;"GOTO" END OF MODIFIED TAYLOR "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";   
        
         "PROCEDURE" DIFFERENCE SCHEME;         
         "BEGIN" "INTEGER" J; "REAL" B;         
             "FOR" J:=M0 "STEP" 1 "UNTIL" M "DO" C[J]:=U[J]; TAUI:=1;         
           NEXT TERM:       
             I:=I+1; DERIVATIVE(I,C); TAUI:=TAUI*TAU;     
             B:=BETA[I]*TAUI;         
             "IF" ETA>0 "AND" I>=P "THEN" LOCAL ERROR CONSTRUCTION(I);        
             "FOR" J:=M0 "STEP" 1 "UNTIL" M "DO" U[J]:=U[J]+B*C[J]; 
             "IF" I<N "THEN" "GOTO" NEXT TERM;  
             T2:=T; "IF" LAST "THEN"  
             "BEGIN" LAST:= "FALSE"; T:= TE "END"         
             "ELSE" T:= T + TAU       
         "END";   
        
         START:= K=0; T0:=T;
         COEFFICIENT; LAST:= "FALSE"; 
       NEXT LEVEL:
         STEPSIZE; K:=K+1; I:=0; DIFFERENCE SCHEME; OUT;  
         "IF" T ^= TE "THEN" "GOTO" NEXT LEVEL  
     "END";       
     END OF MODIFIED TAYLOR:
 "END" MODIFIED TAYLOR;     
         "EOP"    
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 1      
        
        
        
 AUTHORS: P.J. VAN DER HOUWEN AND K.DEKKER.     
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 740416.
        
        
 BRIEF DESCRIPTION:         
        
     EXPONENTIALLY FITTED TAYLOR  SOLVES AN INITIAL VALUE PROBLEM, GIVEN      
     AS A SYSTEM  OF  FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF A      
     ONE-STEP  TAYLOR-METHOD . AUTOMATIC  STEPSIZE  CONTROL IS PROVIDED.      
     IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF  STIFF      
     DIFFERENTIAL EQUATIONS , PROVIDED THAT HIGHER ORDER DERIVATIVES CAN      
     BE EASILY OBTAINED.    
        
        
 KEYWORDS:        
        
     DIFFERENTIAL EQUATIONS,
     INITIAL VALUE PROBLEMS,
     EXPONENTIAL FITTING,   
     STIFF EQUATIONS,       
     THREE-CLUSTER METHOD,  
     ONE-STEP TAYLOR-METHOD.
        
        
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE"  EXPONENTIALLY  FITTED  TAYLOR (T, TE, M0, M, U, SIGMA,      
                            PHI, DIAMETER, DERIVATIVE, I, K, ALFA, NORM,      
                            AETA, RETA, ETA, RHO, HMIN, HSTART, OUTPUT);      
     "INTEGER" M0,M,I,K,NORM;         
     "REAL"  T,TE,SIGMA,PHI,DIAMETER,ALFA,AETA,RETA,ETA,RHO,HMIN,HSTART;      
     "ARRAY" U;   
     "PROCEDURE" DERIVATIVE,OUTPUT;   
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     T:      <VARIABLE>;    
             THE INDEPENDENT VARIABLE T;        
             MAY BE USED IN DERIVATIVE, SIGMA ETC.;       
             ENTRY: THE INITIAL VALUE T0;       
             EXIT : THE FINAL VALUE TE;         
     TE:     <ARITHMETIC EXPRESSION>; 
             THE FINAL VALUE OF T (TE >= T);    
     M0:     <ARITHMETIC EXPRESSION>; 
             INDEX OF THE FIRST EQUATION OF THE SYSTEM TO BE SOLVED;
     M:      <ARITHMETIC EXPRESSION>; 
             INDEX OF THE LAST EQUATION OF THE SYSTEM TO BE SOLVED; 
     U:      <ARRAY IDENTIFIER>;      
             "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 (COMPLEX) POINT  AT WHICH  EXPONENTIAL      
             FITTING  IS DESIRED , FOR EXAMPLE  AN APPROXIMATION  OF THE      
             MODULUS OF THE CENTRE OF THE LEFT HAND CLUSTER;        
     PHI:    <ARITHMETIC EXPRESSION>; 
             THE ARGUMENT  OF THE (COMPLEX) POINT  AT WHICH  EXPONENTIAL      
             FITTING IS DESIRED;      
             PHI SHOULD HAVE A VALUE FROM THE RANGE [PI/2,PI];      
     DIAMETER: <ARITHMETIC EXPRESSION>;         
             THE DIAMETER OF THE LEFT HAND CLUSTER;       
     DERIVATIVE: <PROCEDURE IDENTIFIER>;        
             THE HEADING OF THIS PROCEDURE READS:         
             "PROCEDURE" DERIVATIVE(I,A); "INTEGER" I; "ARRAY" A;   
             I ASSUMES THE VALUES 1,2,3 AND A IS A ONE-DIMENSIONAL ARRAY      
             A[M0:M];       
             WHEN  THIS  PROCEDURE  IS  CALLED , ARRAY  A  CONTAINS  THE      
             COMPONENTS OF THE (I-1)-ST DERIVATIVE OF U  AT THE POINT T;      
             UPON COMPLETION OF DERIVATIVE, ARRAY A SHOULD  CONTAIN  THE      
             COMPONENTS  OF THE I-TH  DERIVATIVE  OF U AT  THE POINT  T;      
     I:      <VARIABLE>;    
             A JENSEN PARAMETER FOR PROCEDURE DERIVATIVE; 
             MAY BE USED IN M0 AND M; 
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 3      
        
        
        
     K:      <VARIABLE>;    
             INDICATES THE NUMBER OF INTEGRATION STEPS PERFORMED;   
             ENTRY: K = 0;  
             EXIT : THE NUMBER OF INTEGRATION STEPS PERFORMED;      
     ALFA:   <ARITHMETIC EXPRESSION>; 
             MAXIMAL GROWTH FACTOR FOR THE INTEGRATION STEP LENGTH; 
     NORM:   <ARITHMETIC EXPRESSION>; 
             IF NORM = 1 DISCREPANCY AND TOLERANCE ARE ESTIMATED IN  THE      
             MAXIMUM NORM, OTHERWISE IN THE EUCLIDIAN NORM;         
     AETA:   <ARITHMETIC EXPRESSION>; 
             DESIRED  ABSOLUTE LOCAL ACCURACY ; AETA SHOULD BE POSITIVE;      
     RETA:   <ARITHMETIC EXPRESSION>; 
             DESIRED  RELATIVE LOCAL ACCURACY ; RETA SHOULD BE POSITIVE;      
     ETA:    <VARIA2LE>;    
             COMPUTED TOLERANCE;      
     RHO:    <VARIABLE>;    
             COMPUTED DISCREPANCY;    
     HMIN:   <ARITHMETIC EXPRESSION>; 
             MINIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
             HOWEVER, A SMALLER STEP  WILL BE TAKEN  IF HMIN EXCEEDS THE      
             STEPSIZE  HSTAB , PRESCRIBED  BY THE  STABILITY  CONDITIONS      
             (SEE REF[2], FORMULA 6.12);        
             IF HSTAB TENDS TO ZERO, THE PROCEDURE TERMINATES;      
     HSTART: <VARIABLE>;    
             ENTRY: THE INTITIAL STEPSIZE ; HOWEVER, IF K = 0  ON ENTRY,      
                    THE VALUE OF HSTART IS NOT TAKEN INTO CONSIDERATION;      
             EXIT:  A SUGGESTION  FOR THE STEPSIZE , IF  THE INTEGRATION      
                    SHOULD BE CONTINUED FOR T>TE;         
             HSTART MAY BE USED IN SUCCESSIVE CALLS OF THE PROCEDURE, IN      
             ORDER TO OBTAIN THE SOLUTION IN SEVERAL POINTS TE1,TE2,ETC;      
     OUTPUT: <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS:         
             "PROCEDURE" OUTPUT;      
             THROUGH  THIS  PROCEDURE  THE VALUES AFTER EACH INTEGRATION      
             STEP  OF  FOR INSTANCE  T, U, ETA AND RHO  ARE  ACCESSIBLE;      
        
        
 DATA AND RESULTS:
        
     FOR FURTHER EXPLANATION OF THE PARAMETERS  SIGMA,PHI,DIAMETER,AETA,      
     RETA,ETA,RHO,M0,M SEE REF[2];    
     FOR RESULTS: SEE EXAMPLE OF USE AND REF[2];
        
        
 PROCEDURES USED: 
        
     INIVEC = CP 31010;     
     DUPVEC = CP 31030;     
     VECVEC = CP 34010;     
     ELMVEC = CP 34020;     
     ZEROIN = CP 34150.     
        
        
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 4      
        
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: CIRCA 40 + 2 * (M - M0).     
        
        
 RUNNING TIME:    
        
     DEPENDS  STRONGLY  ON  THE  DIFFERENTIAL  EQUATION  TO  BE  SOLVED.      
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE REFERENCES.        
        
        
 REFERENCES:      
        
     [1]. P.J. VAN DER HOUWEN.        
          ONE-STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS II,    
          POLYNOMIAL METHODS.         
          TW REPORT 122, (1970) MATHEMATICAL CENTRE.      
        
     [2]. P.J. VAN DER HOUWEN, P.BEENTJES, K.DEKKER AND E.SLAGT.    
          ONE-STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS III,   
          NUMERICAL EXAMPLES.         
          TW REPORT 130, (1971) MATHEMATICAL CENTRE.      
        
        
 EXAMPLE OF USE:  
        
     THE  SOLUTION AT T=EXP(1) AND T=EXP(2) OF THE DIFFERENTIAL EQUATION      
     DU/DT=-EXP(T)*(U-LN(T)) + 1/T WITH INITIAL CONDITION U(.01)=LN(.01)      
     AND  ANALYTICAL SOLUTION U(T) = LN(T), MAY BE OBTAINED  AS FOLLOWS:      
        
     "BEGIN" "INTEGER" I,K; 
         "REAL" T,TE,TE1,TE2,RETA,ETA,RHO,PI,HS,EXPT,LNT,TIME,U0,U1,U2;       
         "REAL" "ARRAY" U[0:0];       
        
         "PROCEDURE" EFT (T,TE,M0,M,U,SIGMA,PHI,DIAMETER,DERIVATIVE,I,K,      
         ALFA,NORM,AETA,RETA,ETA,RHO,HMIN,HSTART,OUTPUT) ; "CODE" 33050;      
        
         "PROCEDURE" DERIVATIVE(I,U); "INTEGER" I; "ARRAY" U;       
         "IF" I=1 "THEN" "BEGIN" EXPT:=EXP(T); LNT:=LN(T); U0:=U[0];
                                 U1:=U[0]:=EXPT*(LNT-U0)+1/T        
                         "END" "ELSE" 
         "IF" I=2 "THEN" U2:=U[0]:=EXPT*(LNT-U0-U1+1/T)-1/T/T       
                  "ELSE" U[0]:=EXPT*(LNT-U0-2*U1-U2+2/T-1/T/T)+2/T/T/T;       
        
         "PROCEDURE" OUT;   
         "IF" T=TE "THEN" OUTPUT(61,"("6ZD,+3ZD.3DB3DB3D")",K,U[0]);
        
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 5      
        
        
        
         "PROCEDURE" OUT1;  
         OUTPUT(61,"("4BD"-D,3Z.3D,/")",RETA,CLOCK-TIME); 
        
         OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
         "("PRINTED BY THIS PROGRAM")",//,      
         "("    THE RESULTS WITH EFT ARE -CONFER REF[2]- :")",/,    
         "("     K        U(TE1)         K        U(TE2)")"         
         "("        RETA   TIME")",/")");       
         PI:=4*ARCTAN(1); TE1:=EXP(1); TE2:=EXP(2);       
         "FOR" RETA:="-1,"-2,"-3,"-4 "DO"       
         "BEGIN" T:=.01; U[0]:=LN(T); K:=0; HS:=0; TIME:=CLOCK;     
             "FOR" TE:=TE1,TE2 "DO"   
             EFT(T,TE,0,0,U,EXP(T),PI,2*EXP(2*T/3),DERIVATIVE,I,K,1.5,2,      
             RETA/10,RETA,ETA,RHO,"-4,HS,OUT); OUT1       
         "END";   
        
         OUTPUT(61,"("//,"("    WITH RELAXED ACCURACY CONDITIONS FOR ")"      
         "("T>3:")",/,"("     K        U(TE1)         K        U(TE2)")"      
         "("        RETA   TIME")",/")");       
         "FOR" RETA:="-1,"-2,"-3,"-4 "DO"       
         "BEGIN" T:=.01; U[0]:=LN(T); K:=0; HS:=0; TIME:=CLOCK;     
             "FOR" TE:=TE1,TE2 "DO"   
             EFT(T,TE,0,0,U,EXP(T),PI,2*EXP(2*T/3),DERIVATIVE,I,K,1.5,2,      
                 RETA/10*("IF" T<3 "THEN" 1 "ELSE" EXP(2*(T-3))),   
                 RETA*("IF" T<3 "THEN" 1 "ELSE" EXP(2*(T-3))),      
                 ETA,RHO,"-4,HS,OUT); OUT1      
         "END"    
     "END"        
        
     THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM    
        
     THE RESULTS WITH EFT ARE -CONFER REF[2]- : 
      K        U(TE1)         K        U(TE2)        RETA   TIME    
      15   +1.003 845 001     42   +2.000 076 417    1"-1   .938    
      22   +1.001 211 286     52   +2.000 066 067    1"-2  1.121    
      36   +1.000 108 738     92   +2.000 020 495    1"-3  1.872    
      56   +1.000 045 271    171   +2.000 000 925    1"-4  3.493    
        
        
     WITH RELAXED ACCURACY CONDITIONS FOR T>3:  
      K        U(TE1)         K        U(TE2)        RETA   TIME    
      15   +1.003 845 001     42   +2.000 076 417    1"-1  1.037    
      22   +1.001 211 286     50   +2.000 049 978    1"-2  1.154    
      36   +1.000 108 738     68   +2.000 023 330    1"-3  1.419    
      56   +1.000 045 271     98   +2.000 065 056    1"-4  2.008    
        
        
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 6      
        
        
        
 SOURCE TEXT(S):  
0"CODE"" 33050;   
 "PROCEDURE" EXPONENTIALLY FITTED TAYLOR(T,TE,M0,M,U,SIGMA,PHI,DIAMETER,      
         DERIVATIVE,I,K,ALFA,NORM,AETA,RETA,ETA,RHO,HMIN,HSTART,OUTPUT);      
 "INTEGER" M0,M,I,K,NORM;   
 "REAL" T,TE,SIGMA,PHI,DIAMETER,ALFA,AETA,RETA,ETA,RHO,HMIN,HSTART; 
 "ARRAY" U;       
 "PROCEDURE" DERIVATIVE,OUTPUT;       
 "BEGIN" "INTEGER" KL;      
     "REAL" Q,EC0,EC1,EC2,H,HI,H0,H1,H2,BETAN,T2,SIGMAL,PHIL;       
     "REAL" "ARRAY" C,RO[M0:M],BETA,BETHA[1:3]; 
     "BOOLEAN" LAST,START;  
     "PROCEDURE" INIVEC(L,U,A,X); "CODE" 31010; 
     "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;   
     "BOOLEAN" "PROCEDURE" ZEROIN(X, Y, FX, EPS); "CODE" 34150;     
        
     "PROCEDURE" COEFFICIENT;         
     "BEGIN" "REAL" B,B1,B2,BB,E,BETA2,BETA3;   
         B:=H*SIGMAL; B1:=B*COS(PHIL); BB:=B*B; 
         "IF" ABS(B)<"-3 "THEN"       
         "BEGIN" BETA2:=.5-BB/24;     
             BETA3:=1/6+B1/12;        
             BETHA[3]:=.5+B1/3        
         "END" "ELSE"       
         "IF" B1<-40 "THEN" 
         "BEGIN" BETA2:=(-2*B1-4*B1*B1/BB+1)/BB;
             BETA3:=(1+2*B1/BB)/BB;   
             BETHA[3]:=1/BB 
         "END" "ELSE"       
         "BEGIN" E:=EXP(B1)/BB; B2:=B*SIN(PHIL);
             BETA2:=(-2*B1-4*B1*B1/BB+1)/BB;    
             BETA3:=(1+2*B1/BB)/BB;   
             "IF" ABS(B2/B)<"-5 "THEN"
             "BEGIN" BETA2:=BETA2-E*(B1-3);     
                 BETA3:=BETA3+E*(B1-2)/B1;      
                 BETHA[3]:=1/BB+E*(B1-1)        
             "END" "ELSE"   
             "BEGIN" BETA2:=BETA2-E*SIN(B2-3*PHIL)/B2*B;  
                 BETA3:=BETA3+E*SIN(B2-2*PHIL)/B2;        
                 BETHA[3]:=1/BB+E*SIN(B2-PHIL)/B2*B;      
             "END"
         "END";   
         BETA[1]:=BETHA[1]:=1;        
         BETA[2]:=BETA2; BETA[3]:=BETA3;        
         BETHA[2]:=1-BB*BETA3; B:=ABS(B);       
         Q:="IF" B<1.5 "THEN" 4-2*B/3 "ELSE" "IF" B<6 "THEN" (30-2*B)/9       
         "ELSE" 2;
     "END";                                                    "COMMENT"      
        
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 7      
                                                                  ; 
        
        
     "REAL" "PROCEDURE" NORMFUNCTION(NORM,W);   
     "INTEGER" NORM; "ARRAY" W;       
     "BEGIN" "INTEGER" J; "REAL" S,X; 
         S:=0;    
         "IF" NORM=1 "THEN" 
         "BEGIN" "FOR" J:=M0 "STEP" 1 "UNTIL" M "DO"      
             "BEGIN" X:=ABS(W[J]); "IF" X>S "THEN" S:=X "END"       
         "END" "ELSE"       
         S:=SQRT(VECVEC(M0,M,0,W,W)); 
         NORMFUNCTION:=S;   
     "END";       
        
     "PROCEDURE" LOCAL ERROR BOUND;   
     ETA:=AETA+RETA * NORMFUNCTION(NORM,U);     
        
     "PROCEDURE" LOCAL ERROR CONSTRUCTION(I); "INTEGER" I;
     "BEGIN" "IF" I=1 "THEN" INIVEC(M0,M,RO,0); 
         "IF" I<4 "THEN" ELMVEC(M0,M,0,RO,C,BETHA[I]*HI); 
         "IF" I=4 "THEN"    
         "BEGIN" ELMVEC(M0,M,0,RO,C,-H);        
             RHO:=NORMFUNCTION(NORM,RO);        
             EC0:=EC1;EC1:=EC2;EC2:=RHO/H**Q;   
         "END"    
     "END";       
        
     "PROCEDURE" STEPSIZE;  
     "BEGIN" "REAL" HACC,HSTAB,HCR,HMAX,A,B,C;  
         "IF" "NOT" START "THEN" LOCAL ERROR BOUND;       
         "IF" START "THEN"  
         "BEGIN" H1:=H2:=HACC:=HSTART;
             EC2:=EC1:=1; KL:=1; START:="FALSE" 
         "END" "ELSE"       
         "IF" KL<3 "THEN"   
         "BEGIN" HACC:=(ETA/RHO)**(1/Q)*H2;     
             "IF" HACC>10*H2 "THEN" HACC:=10*H2 "ELSE" KL:=KL+1     
         "END" "ELSE"       
         "BEGIN" A:=(H0*(EC2-EC1)-H1*(EC1-EC0))/(H2*H0-H1*H1);      
             H:=H2*("IF" ETA<RHO "THEN" (ETA/RHO)**(1/Q) "ELSE" ALFA);        
             "IF" A>0 "THEN"
             "BEGIN" B:=(EC2-EC1-A*(H2-H1))/H1; 
                 C:=EC2-A*H2-B*T2; HACC:=0; HMAX:=H;      
                 "IF" ^ZEROIN(HACC,H,HACC**Q*(A*HACC+B*T+C)-ETA,    
                 "-3*H2) "THEN" HACC:=HMAX      
             "END" "ELSE" HACC:=H;    
             "IF" HACC<.5*H2 "THEN" HACC:=.5*H2;
         "END";   
         "IF" HACC<HMIN "THEN" HACC:=HMIN; H:=HACC;            "COMMENT"      
1SECTION : 5.2.1.1.1.3.B      (AUGUST 1974)                       PAGE 8      
                                                                  ; 
        
        
         "IF" H*SIGMAL>1 "THEN"       
         "BEGIN" A:=ABS(DIAMETER/SIGMAL+"-14)/2; B:=2*ABS(SIN(PHIL));         
             BETAN:=("IF" A>B "THEN" 1/A "ELSE" 1/B)/A;   
             HSTAB:=ABS(BETAN/SIGMAL);
             "IF" HSTAB<"-14*T "THEN" "GOTO" ENDOFEFT;    
             "IF" H>HSTAB "THEN" H:=HSTAB       
         "END";   
         HCR:=H2*H2/H1;     
         "IF" KL>2 "AND" ABS(H-HCR)<"-6*HCR "THEN"        
         H:="IF" H<HCR "THEN" HCR*(1-"-7) "ELSE" HCR*(1+"-7);       
         "IF" T+H>TE "THEN" 
         "BEGIN" LAST:="TRUE"; HSTART:=H; H:=TE-T "END";  
         H0:=H1;H1:=H2;H2:=H;         
     "END";       
        
     "PROCEDURE" DIFFERENCE SCHEME;   
     "BEGIN" HI:=1; SIGMAL:=SIGMA; PHIL:=PHI;   
         STEPSIZE;
         COEFFICIENT;       
         "FOR" I:=1,2,3 "DO"
         "BEGIN" HI:=HI*H;  
             "IF" I>1 "THEN" DERIVATIVE(I,C);   
             LOCALERRORCONSTRUCTION(I);         
             ELMVEC(M0,M,0,U,C,BETA[I]*HI)      
         "END";   
         T2:=T; K:=K+1;     
         "IF" LAST "THEN"   
         "BEGIN" LAST:="FALSE"; T:=TE; START:="TRUE"      
         "END" "ELSE" T:=T+H;         
         DUPVEC(M0,M,0,C,U);
         DERIVATIVE(1,C);   
         LOCALERRORCONSTRUCTION(4);   
         OUTPUT;  
     "END";       
        
     START:="TRUE"; LAST:="FALSE";    
     DUPVEC(M0,M,0,C,U);    
     DERIVATIVE(1,C);       
     "IF" K=0 "THEN"        
     "BEGIN" LOCAL ERROR BOUND; HSTART:=ETA/NORMFUNCTION(NORM,C)    
     "END";       
 NEXT LEVEL:      
     DIFFERENCE SCHEME;     
     "IF" T^=TE "THEN" "GOTO" NEXT LEVEL;       
 ENDOFEFT:        
  "END" EXPONENTIAL FITTED TAYLOR;    
         "EOP"    
1SECTION : 5.2.1.1.2.1        (FEBRUARY 1979)                     PAGE 1
 
 
 
 SECTION 5.2.1.1.2.1 CONTAINS FOUR PROCEDURES FOR INITIAL VALUE PROBLEMS
 FOR SECOND ORDER ORDINARY DIFFERENTIAL EQUATIONS.
 
 A.  RK2 SOLVES AN IVP FOR A SINGLE SECOND ORDER ODE BY MEANS OF A
     5-TH ORDER RUNGE-KUTTA METHOD.
 
 B.  RK2N SOLVES AN IVP FOR A SYSTEM OF SECOND ORDER ODE'S BY MEANS
     OF A 5-TH ORDER RUNGE-KUTTA METHOD
 
 C.  RK3 SOLVES AN IVP FOR A SINGLE SECOND ORDER ODE WITHOUT FIRST
     DERIVATIVE. RK3 IS BASED ON A 5-TH ORDER RUNGE-KUTTA METHOD.
 
 D.  RK3N SOLVES AN IOVP FOR A SYSTEM OF SECOND ORDER ODE'S WITHOUT
     FIRST DERIVATIVE. RK3N IS BASED ON A 5-TH ORDER RUNGE-KUTTA METHOD.
 
1SECTION : 5.2.1.1.2.1.A      (FEBRUARY 1979)                     PAGE 1
 
 
 
 PROCEDURE : RK2.
 
 
 AUTHOR: J.A.ZONNEVELD.
 
 
 CONTRIBUTORS: M.BAKKER AND I.BRINK.
 
 
 INSTITUTE: MATHEMATICAL CENTRE.
 
 
 RECEIVED: 730715.
 
 
 BRIEF DESCRIPTION:
 
     RK2 INTEGRATES THE SCALAR INITIAL VALUE PROBLEM
     (D/DX) (D/DX) Y  = F(X, Y, (D/DX)Y), A<= X <=B OR B <= X <= A,
     Y(A) AND (D/DX) Y(A) PRESCRIBED.
 
 
 KEYWORDS:
 
     INITIAL VALUE PROBLEM,
     SECOND ORDER DIFFERENTIAL EQUATION.
 
 
1SECTION : 5.2.1.1.2.1.A      (FEBRUARY 1979)                     PAGE 2
 
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK2(X, A, B, Y, YA, Z, ZA, FXYZ, E, D, FI);
     "VALUE" B, FI; "REAL" X, A, Y, YA, Z, ZA, FXYZ;
     "BOOLEAN" FI; "ARRAY" E, D;
     "CODE" 33012;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     X:  <VARIABLE>;
         THE INDEPENDENT VARIABLE;
     A:  <ARITHMETIC EXPRESSION>;
         THE INITIAL VALUE OF X;
     B:  <ARITHMETIC EXPRESSION>;
         THE END VALUE OF X, (B <= A IS ALLOWED);
     Y:  <VARIABLE>;
         THE DEPENDENT VARIABLE;
         EXIT :  THE VALUE OF Y(X) AT X = B;
     YA:  <ARITHMETIC EXPRESSION>;
         ENTRY : THE INITIAL VALUE OF Y AT X = A,
     Z:  <VARIABLE>;
         THE DERIVATIVE DY / DX;
         EXIT : THE VALUE OF Z(X) AT X = B;
     ZA:  <ARITHMETIC EXPRESSION>;
         ENTRY : THE INITIAL VALUE OF (D/DX) Y AT X = A;
     FXYZ:  <ARITHMETIC EXPRESSION>;
         THE RIGHT HAND SIDE OF THE DIFFERENTIAL EQUATION;
         FXYZ DEPENDS ON X, Y, Z, GIVING THE VALUE OF (D/DX) (D/DX) Y;
     E:  <ARRAY IDENTIFIER>;
         "ARRAY" E[1 : 4];
         E[1] AND E[3] ARE USED AS RELATIVE , E[2] AND E[4] ARE USED
         AS ABSOLUTE TOLERANCES FOR Y AND DY / DX, RESPECTIVELY;
     D:  <ARRAY IDENTIFIER>;
         "ARRAY" D[1 : 5];
         EXIT:
             ENTIER(D[1] + .5)     =    THE NUMBER OF STEPS SKIPPED,
                    D[2]           =    THE LAST STEP LENGTH USED,
                    D[3]           =    B,
                    D[4]           =    Y(B),
                    D[5]           =    (D/DX) Y, FOR X = B;
     FI:  <BOOLEAN EXPRESSION>;
         IF FI = "TRUE" THEN THE INTEGRATION STARTS AT X=A  WITH A TRIAL
         STEP B - A ; IF FI = "FALSE" THEN THE INTEGRATION IS CONTINUED
         WITH,AS INITIAL CONDITIONS, X = D[3], Y = D[4], Z = D[5], AND
         A, YA AND ZA ARE IGNORED.
 
 
 PROCEDURES USED: NONE.
 
 
1SECTION : 5.2.1.1.2.1.A      (DECEMBER 1979)                     PAGE 3
 
 
 
 METHOD AND PERFORMANCE :
     THE PROCEDURE, WHICH IS PROVIDED WITH STEPLENGTH AND ERROR CONTROL,
     IS BASED ON A 5-TH ORDER RUNGE-KUTTA METHOD.
     A COMPLETE DESCRIPTION IS GIVEN IN [1].
 
 
 REFERENCES:
 
     [1]. J.A.ZONNEVELD.
     AUTOMATIC NUMERICAL INTEGRATION.
     MATH. CENTRE TRACT 8 (1970).
 
 EXAMPLE OF USE:
 
     THE VAN DER POL EQUATION
 
         (D/DX) (D/DX) Y = 10*(1-Y**2)*(DY/DX) - Y, X >= 0,
         Y = 2, DY/DX = 0                , X=0
 
     CAN BE INTEGRATED BY THE PROCEDURE RK2; AT THE POINTS
     X = 9.32386578, 18.86305405, 28.40224162, 37.94142918
     THE DERIVATIVE DY / DX VANISHES; THE PROGRAM WHICH SOLVES THE VAN
     DER POL EQUATION READS AS FOLLOWS (WITH E[I] = "-8, I = 1,...,4):
 
     "BEGIN" "COMMENT" VAN DER POL;
         "PROCEDURE" RK2(X,A,B,Y,YA,Z,ZA,FXYZ,E,D,FI); "CODE" 33012;
         "REAL" X,Y,Z,B; "BOOLEAN" FI; "ARRAY" E[1:4],D[1:5];
         E[1]:=E[2]:=E[3]:=E[4]:="-8;
         "FOR" B:=9.32386578,18.86305405,28.40224162,37.94142918 "DO"
         "BEGIN" FI:= B<10;
             RK2(X,0,B,Y,2,Z,0,10*(1-Y**2)*Z-Y,E,D,FI);
             OUTPUT(61,"("//10B"("X=")"2D.10D,10B"("Y=")"+2D.10D   ,
             10B"("DY/DX =")"+.5D"2D")",X,Y,Z)
         "END"
     "END"
 
     RESULTS:
 
       X=09.3238657800      Y=-02.0142853609      DY/DX=+.00000"00
 
       X=18.8630540500      Y=+02.0142853609      DY/DX=-.00001"00
 
       X=28.4022416200      Y=-02.0142853609      DY/DX=+.00001"00
 
       X=37.9414291800      Y=+02.0142853608      DY/DX=-.00002"00
 
1SECTION : 5.2.1.1.2.1.A      (DECEMBER 1979)                     PAGE 4
 
 
 
 SOURCE TEXT(S):
0"CODE" 33012 ;
    "PROCEDURE" RK2(X, A, B, Y, YA, Z, ZA, FXYZ, E, D, FI);
    "VALUE" B, FI; "REAL" X, A, B, Y, YA, Z, ZA, FXYZ; "BOOLEAN" FI;
    "ARRAY" E, D;
    "BEGIN" "REAL" E1, E2, E3, E4, XL, YL, ZL, H, INT, HMIN, HL,
       ABSH, K0, K1, K2, K3, K4, K5, DISCRY, DISCRZ, TOLY,
       TOLZ, MU, MU1, FHY, FHZ;
       "BOOLEAN" LAST, FIRST, REJECT;
       "IF" FI "THEN"
       "BEGIN" D[3]:= A; D[4]:= YA; D[5]:= ZA "END";
       D[1]:= 0; XL:= D[3]; YL:= D[4]; ZL:= D[5];
       "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]; HL:= INT * E[3] + E[4];
       "IF" HL < HMIN "THEN" HMIN:= HL; E1:= E[1] / INT;
       E2:= E[2] / INT; E3:= E[3] / INT; E4:= E[4] / 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; Z:= ZL; K0:= FXYZ * H;
       X:= XL + H / 4.5;
       Y:= YL + (ZL * 18 + K0 * 2) / 81 * H;
       Z:= ZL + K0 / 4.5  ; K1:= FXYZ * H; X:= XL + H / 3;
       Y:= YL + (ZL * 6 + K0) / 18 * H;
       Z:= ZL + (K0 + K1 * 3) / 12; K2:= FXYZ * H;
       X:= XL + H * .5;
       Y:= YL + (ZL * 8 + K0 + K2) / 16 * H;
       Z:= ZL + (K0 + K2 * 3) / 8; K3:= FXYZ * H;
       X:= XL + H * .8;
       Y:= YL + (ZL * 100 + K0 * 12 + K3 * 28) / 125 * H;      "COMMENT"
1SECTION : 5.2.1.1.2.1.A      (AUGUST 1974)                       PAGE 5
                                                                  ;
 
 
       Z:= ZL + (K0 * 53 - K1 * 135 + K2 * 126 + K3 * 56)
       / 125; K4:= FXYZ * H; X:= "IF" LAST "THEN" B "ELSE" XL + H;
       Y:= YL + (ZL * 336 + K0 * 21 + K2 * 92 + K4 * 55) /
       336 * H;
       Z:= ZL + (K0 * 133 - K1 * 378 + K2 * 276 + K3 * 112
       + K4 * 25) / 168; K5:= FXYZ * H;
       DISCRY:= ABS(( - K0 * 21 + K2 * 108 - K3 * 112 + K4
       * 25) / 56 * H);
       DISCRZ:= ABS(K0 * 21 - K2 * 162 + K3 * 224 - K4 *
       125 + K5 * 42) / 14;
       TOLY:= ABSH * (ABS(ZL) * E1 + E2);
       TOLZ:= ABS(K0) * E3 + ABSH * E4;
       REJECT:= DISCRY > TOLY "OR" DISCRZ > TOLZ;
       FHY:= DISCRY / TOLY; FHZ:= DISCRZ / TOLZ;
       "IF" FHZ > FHY "THEN" FHY:= FHZ;
       MU:= 1 / (1 + FHY) + .45; "IF" REJECT "THEN"
       "BEGIN" "IF" ABSH <= HMIN "THEN"
          "BEGIN" D[1]:= D[1] + 1; Y:= YL; Z:= ZL;
             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";
       FHY:= MU * H / HL + MU - MU1; HL:= H; H:= FHY * H;
    ACC: MU1:= MU;
       Y:= YL + (ZL * 56 + K0 * 7 + K2 * 36 - K4 * 15) / 56
       * HL;
       Z:= ZL + ( - K0 * 63 + K1 * 189 - K2 * 36 - K3 * 112
       + K4 * 50) / 28; K5:= FXYZ * HL;
       Y:= YL + (ZL * 336 + K0 * 35 + K2 * 108 + K4 * 25)
       / 336 * HL;
       Z:= ZL + (K0 * 35 + K2 * 162 + K4 * 125 + K5 * 14)
       / 336;
    NEXT: "IF" B ^= X "THEN"
       "BEGIN" XL:= X; YL:= Y; ZL:= Z; "GOTO" TEST "END";
       "IF" "NOT"LAST "THEN" D[2]:= H; D[3]:= X; D[4]:= Y; D[5]:= Z
    "END" RK2;
         "EOP"
1SECTION : 5.2.1.1.2.1.B      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : RK2N.
        
        
 AUTHOR:J.A.ZONNEVELD.      
        
        
 CONTRIBUTORS: M.BAKKER AND I.BRINK.  
        
        
 INSTITUTE : MATHEMATICAL CENTRE.     
        
        
 RECEIVED: 730715.
        
        
 BRIEF DESCRIPTION: 
        
     RK2N INTEGRATES THE VECTOR INITIAL VALUE PROBLEM     
     (D/DX) (D/DX) Y = F(X, Y, (D/DX) Y), A<= X <= B OR B <= X <= A,
     Y[J] (A)  AND  (D/DX) Y[J] (A) PRESCRIBED FOR J=1,....N.       
        
        
 KEYWORDS :       
        
     INITIAL VALUE PROBLEM, 
     SECOND ORDER DIFFERENTIAL EQUATION.
        
        
1SECTION : 5.2.1.1.2.1.B      (FEBRUARY 1979)                     PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK2N(X,A,B,Y,YA,Z,ZA,FXYZJ,J,E,D,FI,N);  
     "VALUE" B,FI,N;
     "INTEGER" J,N; 
     "REAL" X,A,B,FXYZJ;    
     "BOOLEAN" FI;
     "ARRAY" Y,YA,Z,ZA,E,D; 
     "CODE" 33013;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:    <VARIABLE>;      
           THE INDEPENDENT VARIABLE.  
           UPON COMPLETION OF A CALL OF RK2N,   
           IT IS EQUAL TO B;
     A:    <ARITHMETIC EXPRESSION>;   
           THE STARTING VALUE OF X;   
     B:    <ARITHMETIC EXPRESSION>;   
           A VALUE PARAMETER,GIVING THE END VALUE OF X;   
     Y:    <ARRAY IDENTIFIER>;
           "ARRAY" Y[1:N];  
           THE VECTOR OF DEPENDENT VARIABLES;   
           EXIT: THE VALUE OF Y[J] (B), (J = 1, .. ,N);   
     YA:   <ARRAY IDENTIFIER>;
           "ARRAY" YA[1:N]; 
           ENTRY : THE STARTING VALUES OF Y[J],I.E. THE VALUES AT X=A;
     Z:    <ARRAY IDENTIFIER>;
           "ARRAY" Z[1:N];  
           THE FIRST DERIVATIVES OF THE DEPENDENT VARIABLES;
           EXIT : THE VALUE OF (D/DX)Y[J](B) (J = 1, .. ,N);
     ZA:   <ARRAY IDENTIFIER>;
           "ARRAY" ZA[1:N]; 
           ENTRY : THE STARTING VALUES OF Z[J],I.E. THE VALUES AT X=A;
     FXYZJ:<ARITHMETIC EXPRESSION>;   
           AN EXPRESSION DEPENDING ON X,J,Y[J],Z[J]  (J=1,...,N),   
           GIVING THE VALUE OF (D/DX)(D/DX)Y[J];
     J:    <VARIABLE>;      
           A VARIABLE OF TYPE INTEGER,USED IN THE ACTUAL PARAMETER  
           CORRESPONDING TO FXYZJ,TO DENOTE THE NUMBER OF THE       
           EQUATION REQUIRED (JENSEN'S DEVICE); 
     E:    <ARRAY IDENTIFIER>;
           "ARRAY" E[1:4*N];
           THE ELEMENT E[2*J-1] IS A RELATIVE AND E[2*J] IS AN ABSOLUTE       
           TOLERANCE ASSOCIATED WITH Y[J];      
           E[2*(N+J)-1] IS A RELATIVE AND E[2*(N+J)] IS AN ABSOLUTE 
           TOLERANCE ASSOCIATED WITH Z[J];      
1SECTION : 5.2.1.1.2.1.B      (FEBRUARY 1979)                     PAGE 3      
        
        
        
     D:    <ARRAY IDENTIFIER>;
           "ARRAY" D[1:2*N+3];
           EXIT:  
           ENTIER(D[1]+.5) IS THE NUMBER OF STEPS SKIPPED;
           D[2] IS THE LAST  STEP LENGTH USED;  
           D[3] IS EQUAL TO B;
           D[4],...,D[N+3] ARE EQUAL TO Y[1],...,Y[N] FOR X=B,      
           D[N+4],...,D[2*N+3] ARE EQUAL TO THE DERIVATIVES 
           Z[1],...,Z[N] FOR X=B;     
     FI:   <BOOLEAN EXPRESSION>;      
           IF FI="TRUE" THEN THE INTEGRATION STARTS AT A,WITH A TRIAL 
           STEP B-A;IF FI="FALSE" THEN THE INTEGRATION IS CONTINUED 
           VIZ. WITH INITIAL CONDITIONS:X=D[3],Y[J]=D[J+3],Z[J]=    
           D[N+3+J] AND STEP LENGTH H=D[2]*SIGN(B-D[3]), AND
           A, YA, ZA ARE IGNORED;     
     N:    <ARITHMETIC EXPRESSION>;   
           THE NUMBER OF EQUATIONS.   
        
        
 PROCEDURES USED: NONE.     
        
        
 REQUIRED CENTRAL MEMORY:   
     EIGHT ARRAYS OF ORDER N AND ONE OF ORDER 4 * N ARE USED.       
        
        
 METHOD AND PERFORMANCE :   
     RK2N INTEGRATES (D/DX)(D/DX)Y = F(X,Y,Z) FROM X TO B,WITH, EITHER
     (IF FI = "TRUE") X=A, Y[J]=YA[J], Z[J]=ZA[J], OR  (IF FI="FALSE")
     X = D[3], Y[J]=D[J+3], Z[J]=D[N+J+3], J=1,...,N, USING A 5-TH ORDER      
     RUNGE-KUTTA METHOD.    
     UPON COMPLETION OF A CALL OF RK2N WE HAVE:X=D[3]=B, Y[J]=D[J+3]
     THE VALUE OF THE DEPENDENT VARIABLES FOR X=B, Z[J]=D[N+J+3], THE 
     VALUE OF THE DERIVATIVES OF Y[J] AT X=B, J=1,...,N.  
     RK2N USES AS ITS MINIMAL ABSOLUTE STEP LENGTH
     HMIN=MIN (E[2*J-1]*INT+E[2 *J]) WITH 1<=J<=2*N AND INT=
     ABS(B-("IF" FI "THEN" A "ELSE" D[3])).     
     IF A STEP OF LENGTH ABS(H)<=HMIN IS REJECTED, A STEP SIGN(H)*HMIN
     IS SKIPPED.  A STEP IS REJECTED IF THE ABSOLUTE VALUE OF THE   
     COMPUTED DISCRETIZATION ERROR IS GREATER THAN
     ( ABS(Z[J]) * E[2 * J - 1] + E[2 * J] ) * ABS(H) / INT 
     OR IF THAT TERM IS GREATER THEN (ABS(FXYZJ)*E[2*(J+N)-1
     +E[2*(J+N)])ABS(H)/INT, FOR ANY VALUE OF J ,1<=J<=N (INT=ABS(B-A)).      
     SEE REF[1].  
        
        
1SECTION : 5.2.1.1.2.1.B      (DECEMBER 1975)                     PAGE 4      
        
        
        
 EXAMPLE OF USE:  
        
     THE SECOND ORDER (VECTOR) DIFFERENTIAL EQUATION      
        
       (D/DX)(D/DX)Y[1] = -5*(Y[1] + (D/DX)Y[2]) + Y[2],  
        
       (D/DX)(D/DX)Y[2] = -5*(Y[2] + (D/DX)Y[1]) + Y[1], X>=0,      
        
       Y[1] = (D/DX)Y[2] = 1, Y[2] = (D/DX)Y[1] = 0, X=0  
        
     WITH ANALYTIC SOLUTION 
        
       Y[1] = -EXP(-X)*(EXP(-X)*(EXP(-X)*(EXP(-X)/3+.5)-1)-5/6),    
       Y[2] = -EXP(-X)*(EXP(-X)*(EXP(-X)*(EXP(-X)/3-.5)+1)-5/6)     
        
     CAN BE INTEGRATED BY RK2N FROM 0 TO 5 WITH  1,2,3,4 AS REFERENCE 
     POINTS. THE PROGRAM READS AS FOLLOWS:      
        
     "BEGIN" "REAL" B, X, EXPX; "INTEGER" K; "BOOLEAN" FI;
         "ARRAY" Y,YA,Z,ZA[0:2],E[1:8],D[0:7];  
         "PROCEDURE" RK2N(X,A,B,Y,YA,Z,ZA,FXYZJ,J,E,D,FI,N);"CODE"33013;      
         "FOR"  K:=1,2,3,4,5,6,7,8 "DO" E[K]:="-7;
         YA[1]:=ZA[2]:=1; YA[2]:=ZA[1]:=0; B:=1; AA: FI:=B=1;       
         RK2N(X,0,B,Y,YA,Z,ZA,-5*(Y[K]+Z[K])+("IF"K=1"THEN"Y[2]"ELSE" 
         Y[1]),K,E,D,FI,2); 
         "COMMENT" COMPUTATION  OF THE EXACT VALUES OF Y AND DY/DX; 
         EXPX:=EXP(-X);     
         YA[1]:=-EXPX*(EXPX*(EXPX*(EXPX/3+.5)-1)-5/6);    
         YA[2]:=-EXPX*(EXPX*(EXPX*(EXPX/3-.5)+1)-5/6);    
         ZA[1]:=+EXPX*(EXPX*(EXPX*(EXPX/.75+1.5)-2)-5/6); 
         ZA[2]:=+EXPX*(EXPX*(EXPX*(EXPX/.75-1.5)+2)-5/6); 
         OUTPUT(61,"("/20B"("X=")"D.4D/,
         10B"("Y[1]-YEXACT[1]=")"+.14D ,10B"("Y[2]-YEXACT[2]=")"+.14D4/,      
         10B"("Z[1]-ZEXACT[1]=")"+.14D ,10B"("Z[2]-ZEXACT[2]=")"+.14D 
         5/")",X,Y[1]-YA[1],Y[2]-YA[2],Z[1]-ZA[1],Z[2]-ZA[2]);      
         B:=B+1; "IF" B<5 "THEN" "GO TO" AA     
     "END"
     RESULTS:     
        
           X=1.0000 
      Y[1]-YEXACT[1]=+.00000000002955     Y[2]-YEXACT[2]=+.0000000000567      
      Z[1]-ZEXACT[1]=-.00000000013770     Z[2]-ZEXACT[2]=-.0000000002422      
        
           X=2.0000 
      Y[1]-YEXACT[1]=-.00000000085294     Y[2]-YEXACT[2]=+.0000000001486      
      Z[1]-ZEXACT[1]=+.00000000378800     Z[2]-ZEXACT[2]=-.0000000006509      
        
           X=3.0000 
      Y[1]-YEXACT[1]=-.00000000162707     Y[2]-YEXACT[2]=-.0000000004796      
      Z[1]-ZEXACT[1]=+.00000000803265     Z[2]-ZEXACT[2]=+.0000000019380      
        
           X=4.0000 
      Y[1]-YEXACT[1]=-.00000000117993    Y[2]-YEXACT[2]=-.0000000008505       
      Z[1]-ZEXACT[1]=+.00000000633393     Z[2]-ZEXACT[2]=+.0000000039114      
        
1SECTION : 5.2.1.1.2.1.B      (AUGUST 1974)                       PAGE 5      
        
        
        
 SOURCE TEXT(S):  
0"CODE"" 33013 ;  
    "PROCEDURE" RK2N(X, A, B, Y, YA, Z, ZA, FXYZJ, J, E, D, 
    FI, N); "VALUE" B, FI, N; "INTEGER" J, N; "REAL" X, A, B, FXYZJ;
    "BOOLEAN" FI; "ARRAY" Y, YA, Z, ZA, E, D;   
    "BEGIN" "INTEGER" JJ;   
       "REAL" XL, H, INT, HMIN, HL, ABSH, FHM, DISCRY, DISCRZ,      
       TOLY, TOLZ, MU, MU1, FHY, FHZ; 
       "BOOLEAN" LAST, FIRST, REJECT; 
       "ARRAY" YL, ZL, K0, K1, K2, K3, K4, K5[1:N], EE[1:4 *
       N];
       "IF" FI "THEN"       
       "BEGIN" D[3]:= A;    
          "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"  
          "BEGIN" D[JJ + 3]:= YA[JJ]; D[N + JJ + 3]:= ZA[JJ]
          "END"   
       "END";     
       D[1]:= 0; XL:= D[3]; 
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" YL[JJ]:= D[JJ + 3]; ZL[JJ]:= D[N + JJ + 3] "END";    
       "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];      
       "FOR" JJ:= 2 "STEP" 1 "UNTIL" 2 * N "DO" 
       "BEGIN" HL:= INT * E[2 * JJ - 1] + E[2 * JJ];      
          "IF" HL < HMIN "THEN" HMIN:= HL       
       "END";     
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" 4 * N "DO" EE[JJ]:= E[JJ] / 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:= ABS(H)     
       "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; 
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ]; Z[JJ]:= ZL[JJ] "END";      
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K0[J]:= FXYZJ * H;       
       X:= XL + H / 4.5;    
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ] + (ZL[JJ] * 18 + K0[JJ] * 2) /
          81 * H; Z[JJ]:= ZL[JJ] + K0[JJ] / 4.5;
       "END";                                                  "COMMENT"      
1SECTION : 5.2.1.1.2.1.B      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K1[J]:= FXYZJ * H;       
       X:= XL + H / 3;      
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ] + (ZL[JJ] * 6 + K0[JJ]) / 18 * H;     
          Z[JJ]:= ZL[JJ] + (K0[JJ] + K1[JJ] * 3) / 12     
       "END";     
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K2[J]:= FXYZJ * H;       
       X:= XL + H * .5;     
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ] + (ZL[JJ] * 8 + K0[JJ] + K2[JJ])      
          / 16 * H; 
          Z[JJ]:= ZL[JJ] + (K0[JJ] + K2[JJ] * 3) / 8      
       "END";     
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K3[J]:= FXYZJ * H;       
       X:= XL + H * .8;     
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ] + (ZL[JJ] * 100 + K0[JJ] * 12 +       
          K3[JJ] * 28) / 125 * H;     
          Z[JJ]:= ZL[JJ] + (K0[JJ] * 53 - K1[JJ] * 135 +  
          K2[JJ] * 126 + K3[JJ] * 56) / 125     
       "END";     
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K4[J]:= FXYZJ * H;       
       X:= "IF" LAST "THEN" B "ELSE" XL + H;    
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ] + (ZL[JJ] * 336 + K0[JJ] * 21 +       
          K2[JJ] * 92 + K4[JJ] * 55) / 336 * H; 
          Z[JJ]:= ZL[JJ] + (K0[JJ] * 133 - K1[JJ] * 378 + 
          K2[JJ] * 276 + K3[JJ] * 112 + K4[JJ] * 25) / 168
       "END";     
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K5[J]:= FXYZJ * H;       
       REJECT:= "FALSE"; FHM:= 0;     
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" DISCRY:= ABS(( - K0[JJ] * 21 + K2[JJ] * 108 -
          K3[JJ] * 112 + K4[JJ] * 25) / 56 * H);
          DISCRZ:= ABS(K0[JJ] * 21 - K2[JJ] * 162 + K3[JJ]
          * 224 - K4[JJ] * 125 + K5[JJ] * 42) / 14;       
          TOLY:= ABSH * (ABS(ZL[JJ]) * EE[2 * JJ - 1] +   
          EE[2 * JJ]);      
          TOLZ:= ABS(K0[JJ]) * EE[2 * (JJ + N) - 1] + ABSH
          * EE[2 * (JJ + N)]; 
          REJECT:= DISCRY > TOLY "OR" DISCRZ > TOLZ "OR" REJECT;    
          FHY:= DISCRY / TOLY; FHZ:= DISCRZ / TOLZ;       
          "IF" FHZ > FHY "THEN" FHY:= FHZ;      
          "IF" FHY > FHM "THEN" FHM:= FHY       
       "END";                                                  "COMMENT"      
1SECTION : 5.2.1.1.2.1.B      (AUGUST 1974)                       PAGE 7      
                                                                  ; 
        
        
       MU:= 1 / (1 + FHM) + .45; "IF" REJECT "THEN"       
       "BEGIN" "IF" ABSH <= HMIN "THEN" 
          "BEGIN" D[1]:= D[1] + 1;    
             "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" 
             "BEGIN" Y[JJ]:= YL[JJ]; Z[JJ]:= ZL[JJ] "END";
             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";     
       FHM:= MU * H / HL + MU - MU1; HL:= H; H:= FHM * H; 
    ACC: MU1:= MU;
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ] + (ZL[JJ] * 56 + K0[JJ] * 7 + 
          K2[JJ] * 36 - K4[JJ] * 15) / 56 * HL; 
          Z[JJ]:= ZL[JJ] + ( - K0[JJ] * 63 + K1[JJ] * 189 
          - K2[JJ] * 36 - K3[JJ] * 112 + K4[JJ] * 50) / 28
       "END";     
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K5[J]:= FXYZJ * HL;      
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" Y[JJ]:= YL[JJ] + (ZL[JJ] * 336 + K0[JJ] * 35 +       
          K2[JJ] * 108 + K4[JJ] * 25) / 336 * HL; 
          Z[JJ]:= ZL[JJ] + (K0[JJ] * 35 + K2[JJ] * 162 +  
          K4[JJ] * 125 + K5[JJ] * 14) / 336     
       "END";     
    NEXT: "IF" B ^= X "THEN"
       "BEGIN" XL:= X;      
          "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"  
          "BEGIN" YL[JJ]:= Y[JJ]; ZL[JJ]:= Z[JJ] "END";   
          "GOTO" TEST       
       "END";     
       "IF" "NOT"LAST "THEN" D[2]:= H; D[3]:= X;
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" D[JJ + 3]:= Y[JJ]; D[N + JJ + 3]:= Z[JJ] "END"       
    "END" RK2N;   
         "EOP"    
1SECTION : 5.2.1.1.2.1.C      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : RK3  
        
        
 AUTHOR:J.A.ZONNEVELD.      
        
        
 CONTRIBUTORS: M.BAKKER AND I.BRINK.  
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730715.
        
        
 BRIEF DESCRIPTION: 
        
     RK3 INTEGRATES THE SCALAR INITIAL VALUE PROBLEM      
     (D/DX) (D/DX) Y = F(X,Y)   (WITHOUT THE DERIVATIVE (D/DX) Y IN F),       
     A <= X <= B OR B <= X <= A,  Y(A) AND (D/DX) Y(A) PRESCRIBED.  
        
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEM, 
     SECOND ORDER DIFFERENTIAL EQUATION.
        
        
1SECTION : 5.2.1.1.2.1.C      (DECEMBER 1975)                     PAGE 2      
        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK3(X,A,B,Y,YA,Z,ZA,FXY,E,D,FI); "CODE" 33014;     
     "VALUE" B,FI;
     "REAL" X,A,B,Y,YA,Z,ZA,FXY;      
     "BOOLEAN" FI;
     "ARRAY" E,D; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:    <VARIABLE>;      
           THE INDEPENDENT VARIABLE.  
           UPON COMPLETION OF A CALL OF RK3 ,   
           IT IS EQUAL TO B;
     A:    <ARITHMETIC EXPRESSION>;   
           THE STARTING VALUE OF X;   
     B:    <ARITHMETIC EXPRESSION>;   
           A VALUE PARAMETER, GIVING THE END VALUE OF X;  
           B <= A IS ALLOWED; 
     Y:    <VARIABLE>;      
           THE DEPENDENT VARIABLE;    
           EXIT : THE VALUE OF Y(X) AT X = B;   
     YA:   <ARITHMETIC EXPRESSION>;   
           ENTRY : THE VALUE OF Y AT X=A;       
     Z:    <VARIABLE>;      
           THE DERIVATIVE  DY/DX;     
           EXIT : THE VALUE OF DY/DX AT X = B;  
     ZA:   <ARITHMETIC EXPRESSION>;   
           ENTRY : THE VALUE OF DY/DX AT X=A;   
     FXY:  <ARITHMETIC EXPRESSION>;   
           AN EXPRESSION,DEPENDING ON X AND Y ,GIVING THE VALUE OF  
           (D/DX)(D/DX)Y;   
     E:   <ARRAY IDENTIFIER>; 
           "ARRAY" E[1:4];  
           E[1] AND E[3] ARE USED AS RELATIVE TOLERANCES, 
           E[2] AND E[4] ARE USED AS ABSOLUTE TOLERANCES  
           FOR Y AND DY/DX, RESPECTIVELY;       
     D:    <ARRAY IDENTIFIER>;
           "ARRAY" D[1:5];  
           EXIT:  
           ENTIER(D[1]+.5) IS THE NUMBER OF STEPS SKIPPED;
           D[2] IS THE LAST STEP LENGTH USED;   
           D[3] IS EQUAL TO B;
           D[4] IS EQUAL TO Y(B);     
           D[5] IS EQUAL TO DY/DX FOR X=B;      
     FI:   <BOOLEAN EXPRESSION>;      
           IF FI="TRUE" THEN THE INTEGRATION STARTS AT X=A WITH A TRIAL       
           STEP B-A;IF FI="FALSE" THEN THE INTEGRATION IS CONTINUED 
           VIZ.  WITH THE INITIAL CONDITIONS X=D[3], Y=D[4], Z=D[5] AND       
           STEP LENGTH  H=D[2]*SIGN(B-D[3]); A,YA,ZA ARE IGNORED.   
        
        
1SECTION : 5.2.1.1.2.1.C      (FEBRUARY 1979)                     PAGE 3      
        
        
        
 PROCEDURES USED : NONE.    
        
        
 METHOD AND PERFORMANCE :   
     RK3 INTEGRATES (D/DX)(D/DX)Y = F(X,Y) FROM X TO B,WITH IF FI="TRUE"      
     THEN X=A, Y=YA,DY/DX=ZA ELSE X=D[3], Y=D[4], Z=D[5]. 
     A 5-TH ORDER RUNGE-KUTTA METHOD IS USED.   
     UPON COMPLETION OF A CALL OF RK3 WE HAVE  X=D[3]=B, Y=D[4]=Y[B], 
     Z=D[5], I.E. THE VALUE OF DY/DX FOR X=B.   
     RK3 USES AS ITS MINIMAL ABSOLUTE STEP LENGTH 
     HMIN=MIN (E[2*J-1]*INT+E[2*J]) WITH 1<=J<=2 AND INT= 
     ABS(B-("IF" FI "THEN" A "ELSE" D[3])).     
     IF A STEP OF LENGTH ABS(H)<=HMIN IS REJECTED ,A STEP SIGN(H)*HMIN
     IS SKIPPED.  A STEP IS REJECTED IF THE ABSOLUTE VALUE OF THE LAST
     TERM TAKEN INTO ACCOUNT IS GREATER THEN (ABS(DY/DX)*E[1]+E[2])*
     ABS(H)/INT OR IF THAT TERM IS GREATER THEN (ABS(FXY)*E[3]+E[4])* 
     ABS(H)/INT ( INT = ABS(B - A) ). 
     SEE REF[1].  
        
        
 REFERENCES:      
     [1]J.A.ZONNEVELD.      
        AUTOMATIC NUMERICAL INTEGRATION.
        MATHEMATICAL CENTRE TRACT 8 (1970).     
        
 EXAMPLE OF USE:  
        
     "BEGIN" "COMMENT" SOLUTION OF Y"=X*Y,Y(0)=0,Y'(0)=1; 
     "PROCEDURE" RK3(X,A,B,Y,YA,Z,ZA,FXY,E,D,FI); "CODE" 33014;     
        
     "REAL" "PROCEDURE" YEXACT(X);"VALUE" X;"REAL" X;     
     "BEGIN" "INTEGER" N;"REAL" X3,S,TERM;      
         X3:=X**3;TERM:=X;S:=0;       
         "FOR" N:=3,N+3 "WHILE" ABS(TERM)>"-14 "DO"       
         "BEGIN" S:=S+TERM;TERM:=TERM*X3/N/(N+1)
         "END";   
         YEXACT:=S
     "END";       
        
     "REAL" X,B,Y,Z;"BOOLEAN" FI;"ARRAY" D,E[1:5];
     E[1]:=E[3]:="-8;E[2]:=E[4]:="-12;
     "FOR" B:=.25,.50,.75,1.00 "DO"   
     "BEGIN" FI:=B<.30;     
         RK3(X,0,B,Y,0,Z,1,X*Y,E,D,FI); 
         OUTPUT(61,"("10B"("Y-YEXACT=")".10D"2D,5B"("X=")"Z.2D,     
         5B"("Y=")"2D.10D//")",Y-YEXACT(X),X,Y) 
     "END"
     "END"
1SECTION : 5.2.1.1.2.1.C      (AUGUST 1974)                       PAGE 4      
        
        
        
     DELIVERS:    
        
     Y-YEXACT=0.0000000000     X= .25     Y=00.2503256420 
        
     Y-YEXACT=0.0000000000     X= .50     Y=00.5052238559 
        
     Y-YEXACT=0.0000000000     X= .75     Y=00.7766332813 
        
     Y-YEXACT=0.0000000000     X=1.00     Y=01.0853396481 
        
        
        
 SOURCE TEXT(S):  
0"CODE"" 33014 ;  
    "PROCEDURE" RK3(X, A, B, Y, YA, Z, ZA, FXY, E, D, FI);
    "VALUE" B, FI; "REAL" X, A, B, Y, YA, Z, ZA, FXY; "BOOLEAN" FI; 
    "ARRAY" E, D; 
    "BEGIN" "REAL" E1, E2, E3, E4, XL, YL, ZL, H, INT, HMIN, HL,    
       ABSH, K0, K1, K2, K3, K4, K5, DISCRY, DISCRZ, TOLY,
       TOLZ, MU, MU1, FHY, FHZ;       
       "BOOLEAN" LAST, FIRST, REJECT; 
       "IF" FI "THEN"       
       "BEGIN" D[3]:= A; D[4]:= YA; D[5]:= ZA "END";      
       D[1]:= 0; XL:= D[3]; YL:= D[4]; ZL:= D[5]; 
       "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]; HL:= INT * E[3] + E[4];  
       "IF" HL < HMIN "THEN" HMIN:= HL; E1:= E[1] / INT;  
       E2:= E[2] / INT; E3:= E[3] / INT; E4:= E[4] / INT; 
       FIRST:= REJECT:= "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";                                  "COMMENT"      
1SECTION : 5.2.1.1.2.1.C      (AUGUST 1974)                       PAGE 5      
                                                                  ; 
        
        
    STEP: "IF" REJECT "THEN"
       "BEGIN" X:= XL; Y:= YL; K0:= FXY * H "END" 
       "ELSE" K0:= K5 * H / HL; X:= XL + .276393202250021 * H;      
       Y:= YL + (ZL * .2763932022 50021 + K0 *  
       .038196601125011) * H; K1:= FXY * H;     
       X:= XL + .72360 6797749979 * H;
       Y:= YL + (ZL * .723606797749979 + K1 * .26180      
       3398874989) * H; K2:= FXY * H; X:= XL + H * .5;    
       Y:= YL + (ZL * .5 + K0 * .046875 + K1 *  
       .079824155839840 - K2 * .001699155839840) * H;     
       K4:= FXY * H; X:= "IF" LAST "THEN" B "ELSE" XL + H;
       Y:= YL + (ZL + K0 * .309016994374947 + K2 *
       .190983005625053) * H; K3:= FXY * H;     
       Y:= YL + (ZL + K0 * .083333333333333 + K1 *
       .301502832395825 + K2 * .115163834270842) * H;     
       K5:= FXY * H;
       DISCRY:= ABS(( - K0 * .5 + K1 * 1.809016994374947 +
       K2 * .690983005625053 - K4 * 2) * H);    
       DISCRZ:= ABS((K0 - K3) * 2 - (K1 + K2) * 10 + K4 * 
       16 + K5 * 4); TOLY:= ABSH * (ABS(ZL) * E1 + E2);   
       TOLZ:= ABS(K0) * E3 + ABSH * E4; 
       REJECT:= DISCRY > TOLY "OR" DISCRZ > TOLZ; 
       FHY:= DISCRY / TOLY; FHZ:= DISCRZ / TOLZ;
       "IF" FHZ > FHY "THEN" FHY:= FHZ; 
       MU:= 1 / (1 + FHY) + .45; "IF" REJECT "THEN"       
       "BEGIN" "IF" ABSH <= HMIN "THEN" 
          "BEGIN" D[1]:= D[1] + 1; Y:= YL; Z:= ZL;
             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";     
       FHY:= MU * H / HL + MU - MU1; HL:= H; H:= FHY * H; 
    ACC: MU1:= MU;
       Z:= ZL + (K0 + K3) * .083333333333333 + (K1 + K2) *
       .416666666666667;    
    NEXT: "IF" B ^= X "THEN"
       "BEGIN" XL:= X; YL:= Y; ZL:= Z; "GOTO" TEST "END"; 
       "IF" "NOT"LAST "THEN" D[2]:= H; D[3]:= X; D[4]:= Y; D[5]:= Z 
    "END" RK3;    
         "EOP"    
1SECTION : 5.2.1.1.2.1.D      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 PROCEDURE : RK3N.
        
        
 AUTHOR:J.A.ZONNEVELD.      
        
 CONTRIBUTORS: M.BAKKER AND I.BRINK.  
        
 INSTITUTE:MATHEMATICAL CENTRE.       
        
 RECEIVED: 730715.
        
 BRIEF DESCRIPTION: 
        
     RK3N INTEGRATES THE VECTOR INITIAL VALUE  PROBLEM    
     (D/DX) (D/DX) Y = F(X,Y), A <= X <= B OR B <= X <= A,
     Y[J] (A)  AND  (D/DX) Y[J] (A)  PRESCRIBED.
        
 KEYWORDS:
        
     INITIAL VALUE PROBLEM, 
     SECOND ORDER DIFFERENTIAL EQUATION.
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" RK3N(X,A,B,Y,YA,Z,ZA,FXYJ,J,E,D,FI,N); "CODE" 33015; 
     "VALUE" B,FI,N;
     "INTEGER" J,N; 
     "REAL" X,A,B,FXYJ;     
     "BOOLEAN" FI;
     "ARRAY" Y,YA,Z,ZA,E,D; 
        
1SECTION : 5.2.1.1.2.1.D      (DECEMBER 1975)                     PAGE 2      
        
        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:    <VARIABLE>;      
           THE INDEPENDENT VARIABLE.  
           UPON COMPLETION OF A CALL OF RK3N,   
           IT IS EQUAL TO B;
     A:    <ARITHMETIC EXPRESSION>;   
           THE STARTING VALUE OF X;   
     B:    <ARITHMETIC EXPRESSION>;   
           A VALUE PARAMETER,GIVING THE END VALUE OF X;   
           B <= A IS ALLOWED. 
     Y:    <ARRAY IDENTIFIER>;
           "ARRAY" Y[1:N];  
           THE VECTOR OF DEPENDENT VARIABLES;   
           EXIT : THE VALUE OF Y[J](X) AT X = B;
     YA:   <ARRAY IDENTIFIER>;
           "ARRAY" YA[1:N]; 
           ENTRY : THE STARTING VALUES OF Y[J],I.E. THE VALUES AT X=A;
     Z:    <ARRAY IDENTIFIER>;
           "ARRAY" Z[1:N];  
           THE DERIVATIVES OF THE DEPENDENT VARIABLES, Z[J] = DY[J]/DX;       
           EXIT : THE VALUE OF Z[J](X) AT X = B;
     ZA:   <ARRAY IDENTIFIER>;
           "ARRAY" ZA[1:N]; 
           ENTRY : THE STARTING VALUES OF Z[J],I.E. THE VALUES AT X=A;
     FXYJ: <ARITHMETIC EXPRESSION>;   
           AN EXPRESSION DEPENDING ON X,Y[1],...,Y[N],J,  
           GIVING THE VALUE OF (D/DX)(D/DX)Y[J];
     J:    <VARIABLE>;      
           A VARIABLE OF TYPE INTEGER,USED IN THE ACTUAL PARAMETER  
           CORRESPONDING TO FXYJ,TO DENOTE THE NUMBER OF THE EQUATION 
           REQUIRED (JENSEN'S DEVICE);
     E:    <ARRAY IDENTIFIER>;
           "ARRAY" E[1:4*N];
           THE ELEMENT E[2*J-1] IS A RELATIVE AND E[2*J] IS AN ABSOLUTE       
           TOLERANCE ASSOCIATED WITH Y[J];      
           E[2*(N+J)-1] IS A RELATIVE AND E[2*(N+J)] IS AN ABSOLUTE 
           TOLERANCE ASSOCIATED WITH Z[J];      
     D:    <ARRAY IDENTIFIER>;
           "ARRAY" D[1:2*N+3];
           EXIT:  
           ENTIER(D[1]+.5) IS THE NUMBER OF STEPS SKIPPED;
           D[2] IS THE LAST STEP LENGTH USED;   
           D[3] IS EQUAL TO B;
           D[4],...,D[N+3] ARE EQUAL TO Y[1],...,Y[N] FOR X=B;      
           D[N+4],...,D[2*N+3] ARE EQUAL TO THE DERIVATIVES 
           Z[1],...,Z[N] FOR X=B;     
     FI:   <BOOLEAN EXPRESSION>;      
           IF FI="TRUE" THEN THE INTEGRATION STARTS AT A ,WITH A TRIAL
           STEP B-A;IF FI="FALSE" THEN THE INTEGRATION IS CONTINUED VIZ.      
           WITH THE INITIAL CONDITIONS:X=D[3],Y[J]=D[J+3],Z[J]=D[N+J+3],      
           AND STEP LENGTH H=D[2]*SIGN(B-D[3]); A,YA,ZA ARE IGNORED;
     N:    <ARITHMETIC EXPRESSION>;   
           THE NUMBER OF EQUATIONS.   
        
        
1SECTION : 5.2.1.1.2.1.D      (FEBRUARY 1979)                     PAGE 3      
        
        
        
 PROCEDURES USED: NONE.     
        
 REQUIRED CENTRAL MEMORY:   
     EIGHT ARAYS OF ORDER N AND ONE OF ORDER 4 * N ARE USED.
        
 METHOD AND PERFORMANCE :   
     RK3N INTEGRATES (D/DX)(D/DX)Y=F(X,Y) FROM X TO B,WITH,IF FI="TRUE"       
     THEN X=A, Y[J]=YA[J], Z[J]=ZA[J].IF FI="FALSE" THEN X=D[3],    
     Y[J]=D[J+3], Z[J]=D[N+3+J], USING A 5-TH ORDER RUNGE-KUTTA METHOD.       
     UPON COMPLETION OF A CALL OF RK3N WE HAVE X=D[3]=B, Y[J]=D[J+3]
     THE VALUE OF THE DEPENDENT VARIABLES FOR X=B, Z[J]= D[N+3+J],  
     THE VALUE OF THE DERIVATIVES OF Y[J] AT X=B. 
     RK3N USES AS ITS MINIMAL ABSOLUTE STEP LENGTH:       
     HMIN=MIN (E[2*J-1]*INT+E[2*J]) ,WITH 1<=J<=2*N AND INT=
     ABS(B-("IF" FI "THEN" A "ELSE" D[3])).     
     IF A STEP OF LENGTH ABS(H)<=HMIN IS REJECTED,A STEP SIGN(H)*HMIN IS      
     SKIPPED.     
     A STEP IS REJECTED IF THE ABSOLUTE VALUE OF THE LAST TERM      
     TAKEN INTO ACCOUNT IS GREATER THEN (ABS(Z[J])*E[2*J-1]+E[2*J])*
     ABS(H)/INT OR IF THAT TERM IS GREATER THEN (ABS(FXYJ)*E[2*(J+N)-1]       
     +E[2*(J+N)])*ABS(H)/INT FOR ANY VALUE OF J, 1<=J<=N (INT=ABS(B-A)).      
     SEE REF[1].  
        
 REFERENCES:      
     [1]J.A.ZONNEVELD.      
        AUTOMATIC NUMERICAL INTEGRATION.
        MATHEMATICAL CENTRE TRACT 8 (1970).     
        
 EXAMPLE OF USE:  
     THE SECOND ORDER (VECTOR) DIFFERENTIAL EQUATION      
        
         (D/DX)(D/DX)Y[1] = +Y[2],    
        
         (D/DX)(D/DX)Y[2] = -Y[1], X>=0,
        
         Y[1] = Y[2] = 1,   
        
         (D/DX)Y[1] = (D/DX)Y[2] = 0, X = 0,    
        
     WHOSE EXACT SOLUTION IS GIVEN BY 
        
      Y[1]=COSH(X/SQRT(2))*COS(X/SQRT(2))+SINH(X/SQRT(2))*SIN(X/SQRT(2))      
      Y[2]=COSH(X/SQRT(2))*COS(X/SQRT(2))-SINH(X/SQRT(2))*SIN(X/SQRT(2))      
        
     CAN BE INTEGRATED BY RK3N BECAUSE THE SECOND DERIVATIVE IS NOT 
     EXPRESSED IN THE FIRST. THE PROGRAM READS AS FOLLOWS:
        
1SECTION : 5.2.1.1.2.1.D      (AUGUST 1974)                       PAGE 4      
        
        
        
     "BEGIN" "INTEGER" K,B; "REAL" X; "BOOLEAN" FI;       
         "ARRAY" Y,YA,Z[1:2],E[1:8],D[0:7];     
         "INTEGER" "PROCEDURE" EVEN(N); "VALUE" N; "INTEGER" N;     
         EVEN:= "IF" N//2 = N/2 "THEN" +1 "ELSE" -1;      
         "PROCEDURE" RK3N(X,A,B,Y,YA,Z,ZA,FXYJ,J,E,D,FI,N); "CODE"33015;      
         "PROCEDURE" EXACT(X,Y); "VALUE" X; "REAL" X; "ARRAY" Y;    
         "BEGIN" "INTEGER" I,N; "REAL" X2,TERM; 
             Y[1]:=Y[2]:=0; TERM:=1; X2:= X*X*.5; 
             "FOR" N:=1, N+1 "WHILE" ABS(TERM)>"-14 "DO"  
             "BEGIN" "FOR" I:=1,2 "DO"
                 Y[I]:=Y[I] + TERM*EVEN((I+N-2)//2);      
                 TERM:= TERM*X2  /N/(N*2-1)     
             "END"
         "END";   
         "FOR" K:=1,2,3,4,5,6,7,8 "DO" E[K]:="-7; FI:= "TRUE";      
         Y[1]:=Y[2]:=1; Z[1]:=Z[2]:=0; B:=0; AA: B:= B+1; 
         RK3N(X,0,B,Y,Y,Z,Z,"IF"K=1"THEN"Y[2]"ELSE"-Y[1],K,E,D,FI,2); 
         EXACT(X,YA); OUTPUT(61,"("//10B
         "("ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=")".10D"2D")",  
         ABS(Y[1]-YA[1])+ABS(YA[2]-Y[2]) );     
         FI:="FALSE" ; "IF" B<5 "THEN" "GO TO" AA 
     "END"
     RESULTS:     
     FOR X=1,2,3,4,5 THE FOLLOWING ERRORS ARE NOTICED (E[K]="-7,    
     K=1,...,8):  
           ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000005"00   
           ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000018"00   
           ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000046"00   
           ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000126"00   
           ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000293"00   
        
        
        
 SOURCE TEXT(S):  
0"CODE"" 33015 ;  
    "PROCEDURE" RK3N(X, A, B, Y, YA, Z, ZA, FXYJ, J, E, D,
    FI, N); "VALUE" B, FI, N; "INTEGER" J, N; "REAL" X, A, B, FXYJ; 
    "BOOLEAN" FI; "ARRAY" Y, YA, Z, ZA, E, D;   
    "BEGIN" "INTEGER" JJ;   
       "REAL" XL, H, HMIN, INT, HL, ABSH, FHM, DISCRY, DISCRZ,      
       TOLY, TOLZ, MU, MU1, FHY, FHZ; 
       "BOOLEAN" LAST, FIRST, REJECT; 
       "ARRAY" YL, ZL, K0, K1, K2, K3, K4, K5[1:N], EE[1:4 *
       N];
       "IF" FI "THEN"       
       "BEGIN" D[3]:= A;    
          "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"  
          "BEGIN" D[JJ + 3]:= YA[JJ]; D[N + JJ + 3]:= ZA[JJ]
          "END"   
       "END";                                                  "COMMENT"      
1SECTION : 5.2.1.1.2.1.D      (AUGUST 1974)                       PAGE 5      
                                                                  ; 
        
        
       D[1]:= 0; XL:= D[3]; 
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" YL[JJ]:= D[JJ + 3]; ZL[JJ]:= D[N + JJ + 3] "END";    
       "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];      
       "FOR" JJ:= 2 "STEP" 1 "UNTIL" 2 * N "DO" 
       "BEGIN" HL:= INT * E[2 * JJ - 1] + E[2 * JJ];      
          "IF" HL < HMIN "THEN" HMIN:= HL       
       "END";     
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" 4 * N "DO" EE[JJ]:= E[JJ] / INT; 
       FIRST:= REJECT:= "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: "IF" REJECT "THEN"
       "BEGIN" X:= XL;      
          "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" Y[JJ]:= YL[JJ];      
          "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K0[J]:= FXYJ * H      
       "END"      
       "ELSE"     
       "BEGIN" FHY:= H / HL;
          "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" K0[JJ]:= K5[JJ] * FHY
       "END";     
       X:= XL + .27639 3202250021 * H;
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" Y[JJ]:= YL[JJ] + (ZL[JJ]
       * .276393202250021 + K0[JJ] * .038196601125011) * H; 
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K1[J]:= FXYJ * H;
       X:= XL + .723606797749979 * H; 
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" Y[JJ]:= YL[JJ] + (ZL[JJ]
       * .723606797749979 + K1[JJ] * .261803398874989) * H; 
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K2[J]:= FXYJ * H;
       X:= XL + H * .5;     
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" Y[JJ]:= YL[JJ] + (ZL[JJ]
       * .5 + K0[JJ] * .046875 + K1[JJ] * .079824155839840
       - K2[JJ] * .00169 9155839840) * H;       
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K4[J]:= FXYJ * H;
       X:= "IF" LAST "THEN" B "ELSE" XL + H;    
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" Y[JJ]:= YL[JJ] + (ZL[JJ]
       + K0[JJ] * .309016994374947 + K2[JJ] *   
       .190983005625053) * H;                                  "COMMENT"      
1SECTION : 5.2.1.1.2.1.D      (AUGUST 1974)                       PAGE 6      
                                                                  ; 
        
        
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K3[J]:= FXYJ * H;
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" Y[JJ]:= YL[JJ] + (ZL[JJ]
       + K0[JJ] * .083333333333333 + K1[JJ] * .30150      
       2832395825 + K2[JJ] * .115163834270842) * H;       
       "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" K5[J]:= FXYJ * H;
       REJECT:= "FALSE"; FHM:= 0;     
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" DISCRY:= ABS(( - K0[JJ] * .5 + K1[JJ] *    
          1.809016994374947 + K2[JJ] * .690983005625053 - 
          K4[JJ] * 2) * H); 
          DISCRZ:= ABS((K0[JJ] - K3[JJ]) * 2 - (K1[JJ] +  
          K2[JJ]) * 10 + K4[JJ] * 16 + K5[JJ] * 4);       
          TOLY:= ABSH * (ABS(ZL[JJ]) * EE[2 * JJ - 1] +   
          EE[2 * JJ]);      
          TOLZ:= ABS(K0[JJ]) * EE[2 * (JJ + N) - 1] + ABSH
          * EE[2 * (JJ + N)]; 
          REJECT:= DISCRY > TOLY "OR" DISCRZ > TOLZ "OR" REJECT;    
          FHY:= DISCRY / TOLY; FHZ:= DISCRZ / TOLZ;       
          "IF" FHZ > FHY "THEN" FHY:= FHZ;      
          "IF" FHY > FHM "THEN" FHM:= FHY       
       "END";     
       MU:= 1 / (1 + FHM) + .45; "IF" REJECT "THEN"       
       "BEGIN" "IF" ABSH <= HMIN "THEN" 
          "BEGIN" D[1]:= D[1] + 1;    
             "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" 
             "BEGIN" Y[JJ]:= YL[JJ]; Z[JJ]:= ZL[JJ] "END";
             FIRST:= "TRUE"; "GOTO" NEXT
          "END";  
          H:= MU * H; "GOTO" TEST     
       "END" REJ; 
       "IF" FIRST "THEN"    
       "BEGIN" FIRST:= "FALSE"; HL:= H; H:= MU * H; "GOTO" ACC      
       "END";     
       FHY:= MU * H / HL + MU - MU1; HL:= H; H:= FHY * H; 
    ACC: MU1:= MU;
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO" Z[JJ]:= ZL[JJ] + (K0[JJ]
       + K3[JJ]) * .083333333333333 + (K1[JJ] + K2[JJ]) * 
       .416666666666667;    
    NEXT: "IF" B ^= X "THEN"
       "BEGIN" XL:= X;      
          "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"  
          "BEGIN" YL[JJ]:= Y[JJ]; ZL[JJ]:= Z[JJ] "END";   
          "GOTO" TEST       
       "END";     
       "IF" "NOT"LAST "THEN" D[2]:= H; D[3]:= X;
       "FOR" JJ:= 1 "STEP" 1 "UNTIL" N "DO"     
       "BEGIN" D[JJ + 3]:= Y[JJ]; D[N + JJ + 3]:= Z[JJ] "END"       
    "END" RK3N;   
         "EOP"    
1SECTION : 6.4.1              (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHOR: P.W.HEMKER.
 
 
 CONTRIBUTOR: F.GROEN.
 
 
 INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
 
 
 RECEIVED: 740620.
 
 
 REVISED: 781101 BY N.M.TEMME  AND R.MONTIJN.
 
 
 BRIEF DESCRIPTION:
 
     THIS SECTION CONTAINS THREE PROCEDURES: TAN, ARCSIN, ARCCOS.
 
     TAN COMPUTES THE TANGENT FOR A REAL ARGUMENT X.
     ARCSIN COMPUTES THE ARCSINE FOR A REAL ARGUMENT X.
     ARCCOS COMPUTES THE ARCCOSINE FOR A REAL ARGUMENT X.
 
 KEYWORDS:
 
     TANGENT,
     ARCSINE,
     ARCCOSINE.
1SECTION : 6.4.1              (DECEMBER 1979)                     PAGE 2
 
 
 
 SUBSECTION: TAN.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "REAL" "PROCEDURE" TAN(X); "VALUE" X; "REAL" X;
     "CODE" 35120;
 
     TAN : DELIVERS THE TANGENT OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS:
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF TAN(X).
 
 
 PROCEDURES USED : OVERFLOW  =  CP 30009,
                   GIANT     =  CP 30004.
 
 
 METHOD AND PERFORMANCE :
 
     THE FORMULA TAN(X) = SIN(X) / COS(X) IS USED. IF COS(X) = 0 THEN
     THE VALUE OF GIANT (SEE SECTION 6.2) IS DELIVERED.
 
 EXAMPLE OF USE:
 
     "BEGIN" "REAL""PROCEDURE" TAN(X); "CODE" 35120;
         OUTPUT(61,"("/"("ARCTAN(TAN(1))= ")",+D.14D")",ARCTAN(TAN(1)));
         OUTPUT(61,"("/"("ARCTAN(TAN(0))= ")",+D.14D")",ARCTAN(TAN(0)));
         OUTPUT(61,"("/"("TAN(ARCTAN(0))= ")",+D.14D")",TAN(ARCTAN(0)));
         OUTPUT(61,"("/"("TAN(ARCTAN(1))= ")",+D.14D")",TAN(ARCTAN(1)));
     "END"
 
 DELIVERS :
 
     ARCTAN(TAN(1))= +1.00000000000000
     ARCTAN(TAN(0))= +0.00000000000000
     TAN(ARCTAN(0))= +0.00000000000000
     TAN(ARCTAN(1))= +1.00000000000000
1SECTION : 6.4.1              (DECEMBER 1979)                     PAGE 3
 
 
 
 SUBSECTION : ARCSIN.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "REAL" "PROCEDURE" ARCSIN(X); "VALUE" X; "REAL" X;
     "CODE" 35121;
 
     ARCSIN : DELIVERS THE ARCSINE OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS:
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF ARCSIN(X), ABS(X)<=1.
 
 
 PROCEDURES USED : NONE.
 
 
 METHOD AND PERFORMANCE :
 
     FOR ABS(X) < 0.8 WE USE  THE FORMULA :
         ARCSIN(X) = ARCTAN( X / SQRT ( 1 - X * X )).
     FOR 0.8 <= ABS(X) < 1 WE USE THE FORMULA :
         ARCSIN(X) = SIGN(X) * ( PI/2 - ARCTAN( SQRT( 1/( X * X) - 1))).
     FOR ABS(X) = 1 THE VALUE SIGN(X) * PI/2 IS DELIVERED.
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF ABOUT "-13.
 
 EXAMPLE OF USE :
 
     "BEGIN" "REAL""PROCEDURE" ARCSIN(X); "CODE" 35121;
         OUTPUT(61,"("/"("ARCSIN(SIN(1))= ")",+D.14D")",ARCSIN(SIN(1)));
         OUTPUT(61,"("/"("ARCSIN(SIN(0))= ")",+D.14D")",ARCSIN(SIN(0)));
         OUTPUT(61,"("/"("SIN(ARCSIN(0))= ")",+D.14D")",SIN(ARCSIN(0)));
         OUTPUT(61,"("/"("SIN(ARCSIN(1))= ")",+D.14D")",SIN(ARCSIN(1)));
     "END"
 
     DELIVERS :
 
     ARCSIN(SIN(1))= +0.99999999999990
     ARCSIN(SIN(0))= +0.00000000000000
     SIN(ARCSIN(0))= +0.00000000000000
     SIN(ARCSIN(1))= +1.00000000000000
1SECTION : 6.4.1              (DECEMBER 1979)                     PAGE 4
 
 
 
 SUBSECTION: ARCCOS.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "REAL" "PROCEDURE" ARCCOS(X); "VALUE" X; "REAL" X;
     "CODE" 35122;
 
     ARCCOS : DELIVERS THE ARCCOSINE OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS:
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF ARCCOS(X), ABS(X)<=1.
 
 
 PROCEDURES USED: NONE.
 
 
 METHOD AND PERFORMANCE:
 
     FOR 0 < X < 1 WE USE THE FORMULA:
         ARCCOS(X) = 2 * ARCTAN( SQRT( (1 - X) / (1 + X))).
     FOR -1 < X <= 0 WE USE THE FORMULA:
         ARCCOS(X) = PI - ARCCOS(-X).
     FOR X =  1 THE VALUE 0  IS DELIVERED.
     FOR X = -1 THE VALUE PI IS DELIVERED.
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF "-13.
 
 EXAMPLE OF USE:
 
     "BEGIN" "REAL""PROCEDURE" ARCCOS(X); "CODE" 35122;
         OUTPUT(61,"("/"("ARCCOS(COS(1))= ")",+D.14D")",ARCCOS(COS(1)));
         OUTPUT(61,"("/"("ARCCOS(COS(0))= ")",+D.14D")",ARCCOS(COS(0)));
         OUTPUT(61,"("/"("COS(ARCCOS(0))= ")",+D.14D")",COS(ARCCOS(0)));
         OUTPUT(61,"("/"("COS(ARCCOS(1))= ")",+D.14D")",COS(ARCCOS(1)));
     "END"
 
     DELIVERS :
 
     ARCCOS(COS(1))= +1.00000000000000
     ARCCOS(COS(0))= +0.00000000000000
     COS(ARCCOS(0))= +0.00000000000001
     COS(ARCCOS(1))= +1.00000000000000
 
 
1SECTION : 6.4.1              (DECEMBER 1979)                     PAGE 5
 
 
 
 SOURCE TEXTS:
0"CODE" 35120;
 "REAL" "PROCEDURE" TAN(X); "VALUE" X; "REAL" X;
 "BEGIN" "REAL" U;
     "BOOLEAN" "PROCEDURE" OVERFLOW(X); "CODE" 30009;
     "REAL"    "PROCEDURE" GIANT;       "CODE" 30004;
     U:= SIN(X)/COS(X);
     TAN:= "IF" OVERFLOW(U) "THEN" GIANT "ELSE" U
 "END" TAN;
         "EOP"
 
 "CODE" 35121;
 "REAL" "PROCEDURE" ARCSIN(X); "VALUE" X; "REAL" X;
 "BEGIN" "REAL" U; U:= ABS(X);
     ARCSIN:= "IF" U<0.8 "THEN" ARCTAN(X/SQRT(1-X*X)) "ELSE"
              SIGN(X) * ( "IF" U=1  "THEN" 1.57079632679489 "ELSE"
              ( 1.57079632679489 - ARCTAN(SQRT(1/(X*X)-1))))
 "END" ARCSIN;
         "EOP"
 
 "CODE" 35122;
 "REAL" "PROCEDURE" ARCCOS(X); "VALUE" X; "REAL" X;
 "BEGIN" "REAL" U,V; U:= ABS(X); V:= (1-U)/(1+U);
     V:= "IF" V  =0 "THEN" 0 "ELSE"
         "IF" U+1=1 "THEN" 1.57079632679489 "ELSE"
                           2*ARCTAN(SQRT(V));
     ARCCOS:= "IF" X>0 "THEN" V "ELSE" 3.14159265358979 - V
 "END" ARCCOS;
         "EOP"