share/hompack/fortran/steps.f

   1       SUBROUTINE STEPS(F,NEQN,Y,X,H,EPS,WT,START,HOLD,K,KOLD,CRASH,PHI,
   2      1   P,YP,ALPHA,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI,  FPWA1,FPWA2,FPWA3,
   3      2   FPWA4,FPWA5,IFPWA1,IFPC1,IFPC2,PAR,IPAR)
   4 C
   5 C
   6 C
   7 C   WRITTEN BY L. F. SHAMPINE AND M. K. GORDON
   8 C
   9 C   ABSTRACT
  10 C
  11 C   SUBROUTINE  STEPS  IS NORMALLY USED INDIRECTLY THROUGH SUBROUTINE
  12 C   DEABM .  BECAUSE  DEABM  SUFFICES FOR MOST PROBLEMS AND IS MUCH
  13 C   EASIER TO USE, USING IT SHOULD BE CONSIDERED BEFORE USING  STEPS
  14 C   ALONE.
  15 C
  16 C   SUBROUTINE STEPS INTEGRATES A SYSTEM OF  NEQN  FIRST ORDER ORDINARY
  17 C   DIFFERENTIAL EQUATIONS ONE STEP, NORMALLY FROM X TO X+H, USING A
  18 C   MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE FORMULAS.  LOCAL
  19 C   EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY.
  20 C   THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR
  21 C   PER UNIT STEP IN A GENERALIZED SENSE.  SPECIAL DEVICES ARE INCLUDED
  22 C   TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING
  23 C   TOO MUCH ACCURACY.
  24 C
  25 C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,
  26 C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS, THE INITIAL
  27 C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.
  28 C   FURTHER DETAILS ON USE OF THIS CODE ARE AVAILABLE IN *SOLVING
  29 C   ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
  30 C   BY L. F. SHAMPINE AND M. K. GORDON, SLA-73-1060.
  31 C
  32 C
  33 C   THE PARAMETERS REPRESENT --
  34 C      F -- SUBROUTINE TO EVALUATE DERIVATIVES
  35 C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
  36 C      Y(*) -- SOLUTION VECTOR AT X
  37 C      X -- INDEPENDENT VARIABLE
  38 C      H -- APPROPRIATE STEP SIZE FOR NEXT STEP.  NORMALLY DETERMINED BY
  39 C           CODE
  40 C      EPS -- LOCAL ERROR TOLERANCE
  41 C      WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION
  42 C      START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP,  .FALSE.
  43 C           OTHERWISE
  44 C      HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP
  45 C      K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE)
  46 C      KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP
  47 C      CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN,
  48 C           .FALSE. OTHERWISE.
  49 C      YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT  X  AFTER SUCCESSFUL
  50 C           STEP
  51 C      KSTEPS -- COUNTER ON ATTEMPTED STEPS
  52 C
  53 C   THE VARIABLES X,XOLD,KOLD,KGI AND IVC AND THE ARRAYS Y,PHI,ALPHA,G,
  54 C   W,P,IV AND GI ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE SINTRP.
  55 C   THE ARRAYS FPWA* AND IFPWA1 AND INTEGER CONSTANTS IFPC* ARE
  56 C   WORKING STORAGE PASSED DIRECTLY THROUGH TO  FODE.  THE ARRAYS
  57 C   PAR AND IPAR ARE USER PARAMETERS PASSED THROUGH TO RHOA AND RHOJAC.
  58 C
  59 C   INPUT TO STEPS
  60 C
  61 C      FIRST CALL --
  62 C
  63 C   THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR ALL ARRAYS
  64 C   IN THE CALL LIST, NAMELY
  65 C
  66 C     DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),
  67 C    1  ALPHA(12),W(12),G(13),GI(11),IV(10),   FPWA1(NEQN),
  68 C    2  FPWA2(NEQN-1),FPWA3(NEQN-1,NEQN),FPWA4(NEQN-1),
  69 C    3  FPWA5(NEQN),IFPWA1(NEQN)
  70 C                              --                --    **NOTE**
  71 C
  72 C   THE USER MUST ALSO DECLARE  START  AND  CRASH
  73 C   LOGICAL VARIABLES AND  F  AN EXTERNAL SUBROUTINE, SUPPLY THE
  74 C   SUBROUTINE  F(X,Y,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1,IFPC1,
  75 C                 NEQN-1,IFPC2,PAR,IPAR) TO EVALUATE
  76 C      DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN))
  77 C   AND INITIALIZE ONLY THE FOLLOWING PARAMETERS.
  78 C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
  79 C      Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES
  80 C      X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE
  81 C      H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION
  82 C           AND MAXIMUM SIZE OF STEP.  MUST BE VARIABLE
  83 C      EPS -- LOCAL ERROR TOLERANCE PER STEP.  MUST BE VARIABLE
  84 C      WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION
  85 C      START -- .TRUE.
  86 C      KSTEPS -- SET KSTEPS TO ZERO
  87 C   DEFINE U TO BE THE MACHINE UNIT ROUNDOFF QUANTITY BY CALLING
  88 C   THE FUNCTION ROUTINE  D1MACH,  U = D1MACH(3), OR BY
  89 C   COMPUTING U SO THAT U IS THE SMALLEST POSITIVE NUMBER SUCH
  90 C   THAT 1.0+U .GT. 1.0.
  91 C
  92 C   STEPS  REQUIRES THAT THE L2 NORM OF THE VECTOR WITH COMPONENTS
  93 C   LOCAL ERROR(L)/WT(L)  BE LESS THAN  EPS  FOR A SUCCESSFUL STEP.  THE
  94 C   ARRAY  WT  ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE
  95 C   FOR HIS PROBLEM.  FOR EXAMPLE,
  96 C      WT(L) = 1.0  SPECIFIES ABSOLUTE ERROR,
  97 C            = ABS(Y(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF THE
  98 C                 L-TH COMPONENT OF THE SOLUTION,
  99 C            = ABS(YP(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF
 100 C                 THE L-TH COMPONENT OF THE DERIVATIVE,
 101 C            = MAX(WT(L),ABS(Y(L)))  ERROR RELATIVE TO THE LARGEST
 102 C                 MAGNITUDE OF L-TH COMPONENT OBTAINED SO FAR,
 103 C            = ABS(Y(L))*RELERR/EPS + ABSERR/EPS  SPECIFIES A MIXED
 104 C                 RELATIVE-ABSOLUTE TEST WHERE  RELERR  IS RELATIVE
 105 C                 ERROR,  ABSERR  IS ABSOLUTE ERROR AND  EPS =
 106 C                 MAX(RELERR,ABSERR) .
 107 C
 108 C      SUBSEQUENT CALLS --
 109 C
 110 C   SUBROUTINE  STEPS  IS DESIGNED SO THAT ALL INFORMATION NEEDED TO
 111 C   CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE  H  AND THE ORDER
 112 C   K , IS RETURNED WITH EACH STEP.  WITH THE EXCEPTION OF THE STEP
 113 C   SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS
 114 C   SHOULD BE ALTERED.  THE ARRAY  WT  MUST BE UPDATED AFTER EACH STEP
 115 C   TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE.  NORMALLY THE
 116 C   INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE
 117 C   SOLUTION INTERPOLATED THERE WITH SUBROUTINE  SINTRP .  IF IT IS
 118 C   IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE
 119 C   REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP
 120 C   LARGER THAN THE  H  INPUT.  CHANGING THE DIRECTION OF INTEGRATION,
 121 C   I.E., THE SIGN OF  H , REQUIRES THE USER SET  START = .TRUE. BEFORE
 122 C   CALLING  STEPS  AGAIN.  THIS IS THE ONLY SITUATION IN WHICH  START
 123 C   SHOULD BE ALTERED.
 124 C
 125 C   OUTPUT FROM STEPS
 126 C
 127 C      SUCCESSFUL STEP --
 128 C
 129 C   THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH  START  AND
 130 C   CRASH  SET .FALSE. .  X  REPRESENTS THE INDEPENDENT VARIABLE
 131 C   ADVANCED ONE STEP OF LENGTH  HOLD  FROM ITS VALUE ON INPUT AND  Y
 132 C   THE SOLUTION VECTOR AT THE NEW VALUE OF  X .  ALL OTHER PARAMETERS
 133 C   REPRESENT INFORMATION CORRESPONDING TO THE NEW  X  NEEDED TO
 134 C   CONTINUE THE INTEGRATION.
 135 C
 136 C      UNSUCCESSFUL STEP --
 137 C
 138 C   WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION,
 139 C   THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND  CRASH = .TRUE. .
 140 C   AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE
 141 C   ESTIMATED AND ALL OTHER INFORMATION IS RESTORED AS UPON INPUT
 142 C   BEFORE RETURNING.  TO CONTINUE WITH THE LARGER TOLERANCE, THE USER
 143 C   JUST CALLS THE CODE AGAIN.  A RESTART IS NEITHER REQUIRED NOR
 144 C   DESIRABLE.
 145 C***REFERENCES  SHAMPINE L.F., GORDON M.K., *SOLVING ORDINARY
 146 C                 DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
 147 C                 SLA-73-1060, SANDIA LABORATORIES, 1973.
 148 C
 149       DOUBLE PRECISION ABSH,ALPHA,BETA,D1MACH,EPS,ERK,ERKM1,ERKM2,
 150      1  ERKP1,ERR,FOURU,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,G,GI,GSTR,H,
 151      2  HNEW,HOLD,P,PAR,P5EPS,PHI,PSI,R,REALI,REALNS,RHO,ROUND,SIG,
 152      3  SUM,TAU,TEMP1,TEMP2,TEMP3,TEMP4,TEMP5,TEMP6,TWO,TWOU,V,
 153      4  W,WT,X,XOLD,Y,YP
 154       INTEGER I,IFAIL,IFPC1,IFPC2,IFPWA1,IM1,IPAR,IP1,IQ,IV,IVC,
 155      1  J,JV,K,KGI,KM1,KM2,KNEW,KOLD,KP1,KP2,KPREV,KSTEPS,
 156      2  L,LIMIT1,LIMIT2,NEQN,NS,NSM2,NSP1,NSP2
 157       LOGICAL START,CRASH,PHASE1,NORND
 158 C
 159       DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
 160      1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
 161      2  FPWA1(NEQN),FPWA2(NEQN-1),FPWA3(NEQN-1,NEQN),FPWA4(NEQN-1),
 162      3  FPWA5(NEQN),IFPWA1(NEQN),PAR(1),IPAR(1)
 163       DIMENSION TWO(13),GSTR(13)
 164 C
 165 C   ALL LOCAL VARIABLES ARE SAVED, RATHER THAN PASSED, IN THIS
 166 C   SPECIALIZED VERSION OF STEPS.
 167 C
 168       SAVE
 169 C
 170       EXTERNAL F
 171 C
 172       DATA TWO/2d0,4d0,8d0,16d0,32d0,64d0,128d0,256d0,512d0,1024d0,
 173      1  2048d0,4096d0,8192d0/
 174       DATA GSTR/0.500d0,0.0833d0,0.0417d0,0.0264d0,0.0188d0,0.0143d0,
 175      1  0.0114d0,0.00936d0, 0.00789d0,0.00679d0,0.00592d0,0.00524d0,
 176      2  0.00468d0/
 177 C
 178 C
 179 C       ***     BEGIN BLOCK 0     ***
 180 C   CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE
 181 C   PRECISION.  IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A
 182 C   STARTING STEP SIZE.
 183 C                   ***
 184 C
 185 C   IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE
 186 C
 187 C***FIRST EXECUTABLE STATEMENT
 188       TWOU = 2.0 * D1MACH(4)
 189       FOURU = TWOU + TWOU
 190       CRASH = .TRUE.
 191       IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5
 192       H = SIGN(FOURU*ABS(X),H)
 193       RETURN
 194  5    P5EPS = 0.5*EPS
 195 C
 196 C   IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE
 197 C
 198       ROUND = 0.0
 199       DO 10 L = 1,NEQN
 200  10     ROUND = ROUND + (Y(L)/WT(L))**2
 201       ROUND = TWOU*SQRT(ROUND)
 202       IF(P5EPS .GE. ROUND) GO TO 15
 203       EPS = 2.0*ROUND*(1.0 + FOURU)
 204       RETURN
 205  15   CRASH = .FALSE.
 206       G(1) = 1.0
 207       G(2) = 0.5
 208       SIG(1) = 1.0
 209       IF(.NOT.START) GO TO 99
 210 C
 211 C   INITIALIZE.  COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP
 212 C
 213       CALL F(X,Y,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1,
 214      $       IFPC1,NEQN-1,IFPC2,PAR,IPAR)
 215       IF (IFPC2 .GT. 0) RETURN
 216       SUM = 0.0
 217       DO 20 L = 1,NEQN
 218         PHI(L,1) = YP(L)
 219         PHI(L,2) = 0.0
 220  20     SUM = SUM + (YP(L)/WT(L))**2
 221       SUM = SQRT(SUM)
 222       ABSH = ABS(H)
 223       IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM)
 224       H = SIGN(MAX(ABSH,FOURU*ABS(X)),H)
 225 C
 226 C*      U = D1MACH(3)
 227 C*      BIG = SQRT(D1MACH(2))
 228 C*      CALL HSTART (F,NEQN,X,X+H,Y,YP,WT,1,U,BIG,
 229 C*     1             PHI(1,3),PHI(1,4),PHI(1,5),PHI(1,6),RPAR,IPAR,H)
 230 C
 231       HOLD = 0.0
 232       K = 1
 233       KOLD = 0
 234       KPREV = 0
 235       START = .FALSE.
 236       PHASE1 = .TRUE.
 237       NORND = .TRUE.
 238       IF(P5EPS .GT. 100.0*ROUND) GO TO 99
 239       NORND = .FALSE.
 240       DO 25 L = 1,NEQN
 241  25     PHI(L,15) = 0.0
 242  99   IFAIL = 0
 243 C       ***     END BLOCK 0     ***
 244 C
 245 C       ***     BEGIN BLOCK 1     ***
 246 C   COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP.  AVOID COMPUTING
 247 C   THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.
 248 C                   ***
 249 C
 250  100  KP1 = K+1
 251       KP2 = K+2
 252       KM1 = K-1
 253       KM2 = K-2
 254 C
 255 C   NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT
 256 C   ONE.  WHEN K.LT.NS, NO COEFFICIENTS CHANGE
 257 C
 258       IF(H .NE. HOLD) NS = 0
 259       IF (NS.LE.KOLD) NS = NS+1
 260       NSP1 = NS+1
 261       IF (K .LT. NS) GO TO 199
 262 C
 263 C   COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH
 264 C   ARE CHANGED
 265 C
 266       BETA(NS) = 1.0
 267       REALNS = NS
 268       ALPHA(NS) = 1.0/REALNS
 269       TEMP1 = H*REALNS
 270       SIG(NSP1) = 1.0
 271       IF(K .LT. NSP1) GO TO 110
 272       DO 105 I = NSP1,K
 273         IM1 = I-1
 274         TEMP2 = PSI(IM1)
 275         PSI(IM1) = TEMP1
 276         BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2
 277         TEMP1 = TEMP2 + H
 278         ALPHA(I) = H/TEMP1
 279         REALI = I
 280  105    SIG(I+1) = REALI*ALPHA(I)*SIG(I)
 281  110  PSI(K) = TEMP1
 282 C
 283 C   COMPUTE COEFFICIENTS G(*)
 284 C
 285 C   INITIALIZE V(*) AND SET W(*).
 286 C
 287       IF(NS .GT. 1) GO TO 120
 288       DO 115 IQ = 1,K
 289         TEMP3 = IQ*(IQ+1)
 290         V(IQ) = 1.0/TEMP3
 291  115    W(IQ) = V(IQ)
 292       IVC = 0
 293       KGI = 0
 294       IF (K .EQ. 1) GO TO 140
 295       KGI = 1
 296       GI(1) = W(2)
 297       GO TO 140
 298 C
 299 C   IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)
 300 C
 301  120  IF(K .LE. KPREV) GO TO 130
 302       IF (IVC .EQ. 0) GO TO 122
 303       JV = KP1 - IV(IVC)
 304       IVC = IVC - 1
 305       GO TO 123
 306  122  JV = 1
 307       TEMP4 = K*KP1
 308       V(K) = 1.0/TEMP4
 309       W(K) = V(K)
 310       IF (K .NE. 2) GO TO 123
 311       KGI = 1
 312       GI(1) = W(2)
 313  123  NSM2 = NS-2
 314       IF(NSM2 .LT. JV) GO TO 130
 315       DO 125 J = JV,NSM2
 316         I = K-J
 317         V(I) = V(I) - ALPHA(J+1)*V(I+1)
 318  125    W(I) = V(I)
 319       IF (I .NE. 2) GO TO 130
 320       KGI = NS - 1
 321       GI(KGI) = W(2)
 322 C
 323 C   UPDATE V(*) AND SET W(*)
 324 C
 325  130  LIMIT1 = KP1 - NS
 326       TEMP5 = ALPHA(NS)
 327       DO 135 IQ = 1,LIMIT1
 328         V(IQ) = V(IQ) - TEMP5*V(IQ+1)
 329  135    W(IQ) = V(IQ)
 330       G(NSP1) = W(1)
 331       IF (LIMIT1 .EQ. 1) GO TO 137
 332       KGI = NS
 333       GI(KGI) = W(2)
 334  137  W(LIMIT1+1) = V(LIMIT1+1)
 335       IF (K .GE. KOLD) GO TO 140
 336       IVC = IVC + 1
 337       IV(IVC) = LIMIT1 + 2
 338 C
 339 C   COMPUTE THE G(*) IN THE WORK VECTOR W(*)
 340 C
 341  140  NSP2 = NS + 2
 342       KPREV = K
 343       IF(KP1 .LT. NSP2) GO TO 199
 344       DO 150 I = NSP2,KP1
 345         LIMIT2 = KP2 - I
 346         TEMP6 = ALPHA(I-1)
 347         DO 145 IQ = 1,LIMIT2
 348  145      W(IQ) = W(IQ) - TEMP6*W(IQ+1)
 349  150    G(I) = W(1)
 350  199    CONTINUE
 351 C       ***     END BLOCK 1     ***
 352 C
 353 C       ***     BEGIN BLOCK 2     ***
 354 C   PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED
 355 C   SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,
 356 C   K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.
 357 C                   ***
 358 C
 359 C   INCREMENT COUNTER ON ATTEMPTED STEPS
 360 C
 361       KSTEPS = KSTEPS + 1
 362 C
 363 C   CHANGE PHI TO PHI STAR
 364 C
 365       IF(K .LT. NSP1) GO TO 215
 366       DO 210 I = NSP1,K
 367         TEMP1 = BETA(I)
 368         DO 205 L = 1,NEQN
 369  205      PHI(L,I) = TEMP1*PHI(L,I)
 370  210    CONTINUE
 371 C
 372 C   PREDICT SOLUTION AND DIFFERENCES
 373 C
 374  215  DO 220 L = 1,NEQN
 375         PHI(L,KP2) = PHI(L,KP1)
 376         PHI(L,KP1) = 0.0
 377  220    P(L) = 0.0
 378       DO 230 J = 1,K
 379         I = KP1 - J
 380         IP1 = I+1
 381         TEMP2 = G(I)
 382         DO 225 L = 1,NEQN
 383           P(L) = P(L) + TEMP2*PHI(L,I)
 384  225      PHI(L,I) = PHI(L,I) + PHI(L,IP1)
 385  230    CONTINUE
 386       IF(NORND) GO TO 240
 387       DO 235 L = 1,NEQN
 388         TAU = H*P(L) - PHI(L,15)
 389         P(L) = Y(L) + TAU
 390  235    PHI(L,16) = (P(L) - Y(L)) - TAU
 391       GO TO 250
 392  240  DO 245 L = 1,NEQN
 393  245    P(L) = Y(L) + H*P(L)
 394  250  XOLD = X
 395       X = X + H
 396       ABSH = ABS(H)
 397       CALL F(X,P,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1,
 398      $       IFPC1,NEQN-1,IFPC2,PAR,IPAR)
 399       IF (IFPC2 .GT. 0) RETURN
 400 C
 401 C   ESTIMATE ERRORS AT ORDERS K,K-1,K-2
 402 C
 403       ERKM2 = 0.0
 404       ERKM1 = 0.0
 405       ERK = 0.0
 406       DO 265 L = 1,NEQN
 407         TEMP3 = 1.0/WT(L)
 408         TEMP4 = YP(L) - PHI(L,1)
 409         IF(KM2)265,260,255
 410  255    ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2
 411  260    ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2
 412  265    ERK = ERK + (TEMP4*TEMP3)**2
 413       IF(KM2)280,275,270
 414  270  ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2)
 415  275  ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1)
 416  280  TEMP5 = ABSH*SQRT(ERK)
 417       ERR = TEMP5*(G(K)-G(KP1))
 418       ERK = TEMP5*SIG(KP1)*GSTR(K)
 419       KNEW = K
 420 C
 421 C   TEST IF ORDER SHOULD BE LOWERED
 422 C
 423       IF(KM2)299,290,285
 424  285  IF(MAX(ERKM1,ERKM2) .LE. ERK) KNEW = KM1
 425       GO TO 299
 426  290  IF(ERKM1 .LE. 0.5*ERK) KNEW = KM1
 427 C
 428 C   TEST IF STEP SUCCESSFUL
 429 C
 430  299  IF(ERR .LE. EPS) GO TO 400
 431 C       ***     END BLOCK 2     ***
 432 C
 433 C       ***     BEGIN BLOCK 3     ***
 434 C   THE STEP IS UNSUCCESSFUL.  RESTORE  X, PHI(*,*), PSI(*) .
 435 C   IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE.  IF STEP FAILS MORE
 436 C   THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE.  DOUBLE ERROR
 437 C   TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE
 438 C   PRECISION.
 439 C                   ***
 440 C
 441 C   RESTORE X, PHI(*,*) AND PSI(*)
 442 C
 443       PHASE1 = .FALSE.
 444       X = XOLD
 445       DO 310 I = 1,K
 446         TEMP1 = 1.0/BETA(I)
 447         IP1 = I+1
 448         DO 305 L = 1,NEQN
 449  305      PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1))
 450  310    CONTINUE
 451       IF(K .LT. 2) GO TO 320
 452       DO 315 I = 2,K
 453  315    PSI(I-1) = PSI(I) - H
 454 C
 455 C   ON THIRD FAILURE, SET ORDER TO ONE.  THEREAFTER, USE OPTIMAL STEP
 456 C   SIZE
 457 C
 458  320  IFAIL = IFAIL + 1
 459       TEMP2 = 0.5
 460       IF(IFAIL - 3) 335,330,325
 461  325  IF(P5EPS .LT. 0.25*ERK) TEMP2 = SQRT(P5EPS/ERK)
 462  330  KNEW = 1
 463  335  H = TEMP2*H
 464       K = KNEW
 465       NS = 0
 466       IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340
 467       CRASH = .TRUE.
 468       H = SIGN(FOURU*ABS(X),H)
 469       EPS = EPS + EPS
 470       RETURN
 471  340  GO TO 100
 472 C       ***     END BLOCK 3     ***
 473 C
 474 C       ***     BEGIN BLOCK 4     ***
 475 C   THE STEP IS SUCCESSFUL.  CORRECT THE PREDICTED SOLUTION, EVALUATE
 476 C   THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE
 477 C   DIFFERENCES.  DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.
 478 C                   ***
 479  400  KOLD = K
 480       HOLD = H
 481 C
 482 C   CORRECT AND EVALUATE
 483 C
 484       TEMP1 = H*G(KP1)
 485       IF(NORND) GO TO 410
 486       DO 405 L = 1,NEQN
 487         TEMP3 = Y(L)
 488         RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16)
 489         Y(L) = P(L) + RHO
 490         PHI(L,15) = (Y(L) - P(L)) - RHO
 491  405    P(L) = TEMP3
 492       GO TO 420
 493  410  DO 415 L = 1,NEQN
 494         TEMP3 = Y(L)
 495         Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1))
 496  415    P(L) = TEMP3
 497  420  CALL F(X,Y,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1,
 498      $       IFPC1,NEQN-1,IFPC2,PAR,IPAR)
 499       IF (IFPC2 .GT. 0) RETURN
 500 C
 501 C   UPDATE DIFFERENCES FOR NEXT STEP
 502 C
 503       DO 425 L = 1,NEQN
 504         PHI(L,KP1) = YP(L) - PHI(L,1)
 505  425    PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2)
 506       DO 435 I = 1,K
 507         DO 430 L = 1,NEQN
 508  430      PHI(L,I) = PHI(L,I) + PHI(L,KP1)
 509  435    CONTINUE
 510 C
 511 C   ESTIMATE ERROR AT ORDER K+1 UNLESS:
 512 C     IN FIRST PHASE WHEN ALWAYS RAISE ORDER,
 513 C     ALREADY DECIDED TO LOWER ORDER,
 514 C     STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE
 515 C
 516       ERKP1 = 0.0
 517       IF(KNEW .EQ. KM1  .OR.  K .EQ. 12) PHASE1 = .FALSE.
 518       IF(PHASE1) GO TO 450
 519       IF(KNEW .EQ. KM1) GO TO 455
 520       IF(KP1 .GT. NS) GO TO 460
 521       DO 440 L = 1,NEQN
 522  440    ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2
 523       ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1)
 524 C
 525 C   USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER
 526 C   FOR NEXT STEP
 527 C
 528       IF(K .GT. 1) GO TO 445
 529       IF(ERKP1 .GE. 0.5*ERK) GO TO 460
 530       GO TO 450
 531  445  IF(ERKM1 .LE. MIN(ERK,ERKP1)) GO TO 455
 532       IF(ERKP1 .GE. ERK  .OR.  K .EQ. 12) GO TO 460
 533 C
 534 C   HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE
 535 C   BEEN LOWERED IN BLOCK 2.  THUS ORDER IS TO BE RAISED
 536 C
 537 C   RAISE ORDER
 538 C
 539  450  K = KP1
 540       ERK = ERKP1
 541       GO TO 460
 542 C
 543 C   LOWER ORDER
 544 C
 545  455  K = KM1
 546       ERK = ERKM1
 547 C
 548 C   WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP
 549 C
 550  460  HNEW = H + H
 551       IF(PHASE1) GO TO 465
 552       IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465
 553       HNEW = H
 554       IF(P5EPS .GE. ERK) GO TO 465
 555       TEMP2 = K+1
 556       R = (P5EPS/ERK)**(1.0/TEMP2)
 557       HNEW = ABSH*MAX(0.5D0,MIN(0.9D0,R))
 558       HNEW = SIGN(MAX(HNEW,FOURU*ABS(X)),H)
 559  465  H = HNEW
 560       RETURN
 561 C       ***     END BLOCK 4     ***
 562       END