SUBROUTINE STEPS(F,NEQN,Y,X,H,EPS,WT,START,HOLD,K,KOLD,CRASH,PHI, 1 P,YP,ALPHA,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI, FPWA1,FPWA2,FPWA3, 2 FPWA4,FPWA5,IFPWA1,IFPC1,IFPC2,PAR,IPAR) C C C C WRITTEN BY L. F. SHAMPINE AND M. K. GORDON C C ABSTRACT C C SUBROUTINE STEPS IS NORMALLY USED INDIRECTLY THROUGH SUBROUTINE C DEABM . BECAUSE DEABM SUFFICES FOR MOST PROBLEMS AND IS MUCH C EASIER TO USE, USING IT SHOULD BE CONSIDERED BEFORE USING STEPS C ALONE. C C SUBROUTINE STEPS INTEGRATES A SYSTEM OF NEQN FIRST ORDER ORDINARY C DIFFERENTIAL EQUATIONS ONE STEP, NORMALLY FROM X TO X+H, USING A C MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE FORMULAS. LOCAL C EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY. C THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR C PER UNIT STEP IN A GENERALIZED SENSE. SPECIAL DEVICES ARE INCLUDED C TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING C TOO MUCH ACCURACY. C C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS, THE INITIAL C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. C FURTHER DETAILS ON USE OF THIS CODE ARE AVAILABLE IN *SOLVING C ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*, C BY L. F. SHAMPINE AND M. K. GORDON, SLA-73-1060. C C C THE PARAMETERS REPRESENT -- C F -- SUBROUTINE TO EVALUATE DERIVATIVES C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C Y(*) -- SOLUTION VECTOR AT X C X -- INDEPENDENT VARIABLE C H -- APPROPRIATE STEP SIZE FOR NEXT STEP. NORMALLY DETERMINED BY C CODE C EPS -- LOCAL ERROR TOLERANCE C WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION C START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP, .FALSE. C OTHERWISE C HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP C K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE) C KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP C CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN, C .FALSE. OTHERWISE. C YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT X AFTER SUCCESSFUL C STEP C KSTEPS -- COUNTER ON ATTEMPTED STEPS C C THE VARIABLES X,XOLD,KOLD,KGI AND IVC AND THE ARRAYS Y,PHI,ALPHA,G, C W,P,IV AND GI ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE SINTRP. C THE ARRAYS FPWA* AND IFPWA1 AND INTEGER CONSTANTS IFPC* ARE C WORKING STORAGE PASSED DIRECTLY THROUGH TO FODE. THE ARRAYS C PAR AND IPAR ARE USER PARAMETERS PASSED THROUGH TO RHOA AND RHOJAC. C C INPUT TO STEPS C C FIRST CALL -- C C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR ALL ARRAYS C IN THE CALL LIST, NAMELY C C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN), C 1 ALPHA(12),W(12),G(13),GI(11),IV(10), FPWA1(NEQN), C 2 FPWA2(NEQN-1),FPWA3(NEQN-1,NEQN),FPWA4(NEQN-1), C 3 FPWA5(NEQN),IFPWA1(NEQN) C -- -- **NOTE** C C THE USER MUST ALSO DECLARE START AND CRASH C LOGICAL VARIABLES AND F AN EXTERNAL SUBROUTINE, SUPPLY THE C SUBROUTINE F(X,Y,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1,IFPC1, C NEQN-1,IFPC2,PAR,IPAR) TO EVALUATE C DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN)) C AND INITIALIZE ONLY THE FOLLOWING PARAMETERS. C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES C X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE C H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION C AND MAXIMUM SIZE OF STEP. MUST BE VARIABLE C EPS -- LOCAL ERROR TOLERANCE PER STEP. MUST BE VARIABLE C WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION C START -- .TRUE. C KSTEPS -- SET KSTEPS TO ZERO C DEFINE U TO BE THE MACHINE UNIT ROUNDOFF QUANTITY BY CALLING C THE FUNCTION ROUTINE D1MACH, U = D1MACH(3), OR BY C COMPUTING U SO THAT U IS THE SMALLEST POSITIVE NUMBER SUCH C THAT 1.0+U .GT. 1.0. C C STEPS REQUIRES THAT THE L2 NORM OF THE VECTOR WITH COMPONENTS C LOCAL ERROR(L)/WT(L) BE LESS THAN EPS FOR A SUCCESSFUL STEP. THE C ARRAY WT ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE C FOR HIS PROBLEM. FOR EXAMPLE, C WT(L) = 1.0 SPECIFIES ABSOLUTE ERROR, C = ABS(Y(L)) ERROR RELATIVE TO THE MOST RECENT VALUE OF THE C L-TH COMPONENT OF THE SOLUTION, C = ABS(YP(L)) ERROR RELATIVE TO THE MOST RECENT VALUE OF C THE L-TH COMPONENT OF THE DERIVATIVE, C = MAX(WT(L),ABS(Y(L))) ERROR RELATIVE TO THE LARGEST C MAGNITUDE OF L-TH COMPONENT OBTAINED SO FAR, C = ABS(Y(L))*RELERR/EPS + ABSERR/EPS SPECIFIES A MIXED C RELATIVE-ABSOLUTE TEST WHERE RELERR IS RELATIVE C ERROR, ABSERR IS ABSOLUTE ERROR AND EPS = C MAX(RELERR,ABSERR) . C C SUBSEQUENT CALLS -- C C SUBROUTINE STEPS IS DESIGNED SO THAT ALL INFORMATION NEEDED TO C CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE H AND THE ORDER C K , IS RETURNED WITH EACH STEP. WITH THE EXCEPTION OF THE STEP C SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS C SHOULD BE ALTERED. THE ARRAY WT MUST BE UPDATED AFTER EACH STEP C TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE. NORMALLY THE C INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE C SOLUTION INTERPOLATED THERE WITH SUBROUTINE SINTRP . IF IT IS C IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE C REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP C LARGER THAN THE H INPUT. CHANGING THE DIRECTION OF INTEGRATION, C I.E., THE SIGN OF H , REQUIRES THE USER SET START = .TRUE. BEFORE C CALLING STEPS AGAIN. THIS IS THE ONLY SITUATION IN WHICH START C SHOULD BE ALTERED. C C OUTPUT FROM STEPS C C SUCCESSFUL STEP -- C C THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH START AND C CRASH SET .FALSE. . X REPRESENTS THE INDEPENDENT VARIABLE C ADVANCED ONE STEP OF LENGTH HOLD FROM ITS VALUE ON INPUT AND Y C THE SOLUTION VECTOR AT THE NEW VALUE OF X . ALL OTHER PARAMETERS C REPRESENT INFORMATION CORRESPONDING TO THE NEW X NEEDED TO C CONTINUE THE INTEGRATION. C C UNSUCCESSFUL STEP -- C C WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION, C THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND CRASH = .TRUE. . C AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE C ESTIMATED AND ALL OTHER INFORMATION IS RESTORED AS UPON INPUT C BEFORE RETURNING. TO CONTINUE WITH THE LARGER TOLERANCE, THE USER C JUST CALLS THE CODE AGAIN. A RESTART IS NEITHER REQUIRED NOR C DESIRABLE. C***REFERENCES SHAMPINE L.F., GORDON M.K., *SOLVING ORDINARY C DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*, C SLA-73-1060, SANDIA LABORATORIES, 1973. C DOUBLE PRECISION ABSH,ALPHA,BETA,D1MACH,EPS,ERK,ERKM1,ERKM2, 1 ERKP1,ERR,FOURU,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,G,GI,GSTR,H, 2 HNEW,HOLD,P,PAR,P5EPS,PHI,PSI,R,REALI,REALNS,RHO,ROUND,SIG, 3 SUM,TAU,TEMP1,TEMP2,TEMP3,TEMP4,TEMP5,TEMP6,TWO,TWOU,V, 4 W,WT,X,XOLD,Y,YP INTEGER I,IFAIL,IFPC1,IFPC2,IFPWA1,IM1,IPAR,IP1,IQ,IV,IVC, 1 J,JV,K,KGI,KM1,KM2,KNEW,KOLD,KP1,KP2,KPREV,KSTEPS, 2 L,LIMIT1,LIMIT2,NEQN,NS,NSM2,NSP1,NSP2 LOGICAL START,CRASH,PHASE1,NORND C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12), 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10), 2 FPWA1(NEQN),FPWA2(NEQN-1),FPWA3(NEQN-1,NEQN),FPWA4(NEQN-1), 3 FPWA5(NEQN),IFPWA1(NEQN),PAR(1),IPAR(1) DIMENSION TWO(13),GSTR(13) C C ALL LOCAL VARIABLES ARE SAVED, RATHER THAN PASSED, IN THIS C SPECIALIZED VERSION OF STEPS. C SAVE C EXTERNAL F C DATA TWO/2.0,4.0,8.0,16.0,32.0,64.0,128.0,256.0,512.0,1024.0, 1 2048.0,4096.0,8192.0/ DATA GSTR/0.500,0.0833,0.0417,0.0264,0.0188,0.0143,0.0114,0.00936, 1 0.00789,0.00679,0.00592,0.00524,0.00468/ C C C *** BEGIN BLOCK 0 *** C CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE C PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A C STARTING STEP SIZE. C *** C C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE C C***FIRST EXECUTABLE STATEMENT TWOU = 2.0 * D1MACH(4) FOURU = TWOU + TWOU CRASH = .TRUE. IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5 H = SIGN(FOURU*ABS(X),H) RETURN 5 P5EPS = 0.5*EPS C C IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE C ROUND = 0.0 DO 10 L = 1,NEQN 10 ROUND = ROUND + (Y(L)/WT(L))**2 ROUND = TWOU*SQRT(ROUND) IF(P5EPS .GE. ROUND) GO TO 15 EPS = 2.0*ROUND*(1.0 + FOURU) RETURN 15 CRASH = .FALSE. G(1) = 1.0 G(2) = 0.5 SIG(1) = 1.0 IF(.NOT.START) GO TO 99 C C INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP C CALL F(X,Y,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1, $ IFPC1,NEQN-1,IFPC2,PAR,IPAR) IF (IFPC2 .GT. 0) RETURN SUM = 0.0 DO 20 L = 1,NEQN PHI(L,1) = YP(L) PHI(L,2) = 0.0 20 SUM = SUM + (YP(L)/WT(L))**2 SUM = SQRT(SUM) ABSH = ABS(H) IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM) H = SIGN(MAX(ABSH,FOURU*ABS(X)),H) C C* U = D1MACH(3) C* BIG = SQRT(D1MACH(2)) C* CALL HSTART (F,NEQN,X,X+H,Y,YP,WT,1,U,BIG, C* 1 PHI(1,3),PHI(1,4),PHI(1,5),PHI(1,6),RPAR,IPAR,H) C HOLD = 0.0 K = 1 KOLD = 0 KPREV = 0 START = .FALSE. PHASE1 = .TRUE. NORND = .TRUE. IF(P5EPS .GT. 100.0*ROUND) GO TO 99 NORND = .FALSE. DO 25 L = 1,NEQN 25 PHI(L,15) = 0.0 99 IFAIL = 0 C *** END BLOCK 0 *** C C *** BEGIN BLOCK 1 *** C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID COMPUTING C THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED. C *** C 100 KP1 = K+1 KP2 = K+2 KM1 = K-1 KM2 = K-2 C C NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT C ONE. WHEN K.LT.NS, NO COEFFICIENTS CHANGE C IF(H .NE. HOLD) NS = 0 IF (NS.LE.KOLD) NS = NS+1 NSP1 = NS+1 IF (K .LT. NS) GO TO 199 C C COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH C ARE CHANGED C BETA(NS) = 1.0 REALNS = NS ALPHA(NS) = 1.0/REALNS TEMP1 = H*REALNS SIG(NSP1) = 1.0 IF(K .LT. NSP1) GO TO 110 DO 105 I = NSP1,K IM1 = I-1 TEMP2 = PSI(IM1) PSI(IM1) = TEMP1 BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2 TEMP1 = TEMP2 + H ALPHA(I) = H/TEMP1 REALI = I 105 SIG(I+1) = REALI*ALPHA(I)*SIG(I) 110 PSI(K) = TEMP1 C C COMPUTE COEFFICIENTS G(*) C C INITIALIZE V(*) AND SET W(*). C IF(NS .GT. 1) GO TO 120 DO 115 IQ = 1,K TEMP3 = IQ*(IQ+1) V(IQ) = 1.0/TEMP3 115 W(IQ) = V(IQ) IVC = 0 KGI = 0 IF (K .EQ. 1) GO TO 140 KGI = 1 GI(1) = W(2) GO TO 140 C C IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*) C 120 IF(K .LE. KPREV) GO TO 130 IF (IVC .EQ. 0) GO TO 122 JV = KP1 - IV(IVC) IVC = IVC - 1 GO TO 123 122 JV = 1 TEMP4 = K*KP1 V(K) = 1.0/TEMP4 W(K) = V(K) IF (K .NE. 2) GO TO 123 KGI = 1 GI(1) = W(2) 123 NSM2 = NS-2 IF(NSM2 .LT. JV) GO TO 130 DO 125 J = JV,NSM2 I = K-J V(I) = V(I) - ALPHA(J+1)*V(I+1) 125 W(I) = V(I) IF (I .NE. 2) GO TO 130 KGI = NS - 1 GI(KGI) = W(2) C C UPDATE V(*) AND SET W(*) C 130 LIMIT1 = KP1 - NS TEMP5 = ALPHA(NS) DO 135 IQ = 1,LIMIT1 V(IQ) = V(IQ) - TEMP5*V(IQ+1) 135 W(IQ) = V(IQ) G(NSP1) = W(1) IF (LIMIT1 .EQ. 1) GO TO 137 KGI = NS GI(KGI) = W(2) 137 W(LIMIT1+1) = V(LIMIT1+1) IF (K .GE. KOLD) GO TO 140 IVC = IVC + 1 IV(IVC) = LIMIT1 + 2 C C COMPUTE THE G(*) IN THE WORK VECTOR W(*) C 140 NSP2 = NS + 2 KPREV = K IF(KP1 .LT. NSP2) GO TO 199 DO 150 I = NSP2,KP1 LIMIT2 = KP2 - I TEMP6 = ALPHA(I-1) DO 145 IQ = 1,LIMIT2 145 W(IQ) = W(IQ) - TEMP6*W(IQ+1) 150 G(I) = W(1) 199 CONTINUE C *** END BLOCK 1 *** C C *** BEGIN BLOCK 2 *** C PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED C SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K, C K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED. C *** C C INCREMENT COUNTER ON ATTEMPTED STEPS C KSTEPS = KSTEPS + 1 C C CHANGE PHI TO PHI STAR C IF(K .LT. NSP1) GO TO 215 DO 210 I = NSP1,K TEMP1 = BETA(I) DO 205 L = 1,NEQN 205 PHI(L,I) = TEMP1*PHI(L,I) 210 CONTINUE C C PREDICT SOLUTION AND DIFFERENCES C 215 DO 220 L = 1,NEQN PHI(L,KP2) = PHI(L,KP1) PHI(L,KP1) = 0.0 220 P(L) = 0.0 DO 230 J = 1,K I = KP1 - J IP1 = I+1 TEMP2 = G(I) DO 225 L = 1,NEQN P(L) = P(L) + TEMP2*PHI(L,I) 225 PHI(L,I) = PHI(L,I) + PHI(L,IP1) 230 CONTINUE IF(NORND) GO TO 240 DO 235 L = 1,NEQN TAU = H*P(L) - PHI(L,15) P(L) = Y(L) + TAU 235 PHI(L,16) = (P(L) - Y(L)) - TAU GO TO 250 240 DO 245 L = 1,NEQN 245 P(L) = Y(L) + H*P(L) 250 XOLD = X X = X + H ABSH = ABS(H) CALL F(X,P,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1, $ IFPC1,NEQN-1,IFPC2,PAR,IPAR) IF (IFPC2 .GT. 0) RETURN C C ESTIMATE ERRORS AT ORDERS K,K-1,K-2 C ERKM2 = 0.0 ERKM1 = 0.0 ERK = 0.0 DO 265 L = 1,NEQN TEMP3 = 1.0/WT(L) TEMP4 = YP(L) - PHI(L,1) IF(KM2)265,260,255 255 ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2 260 ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2 265 ERK = ERK + (TEMP4*TEMP3)**2 IF(KM2)280,275,270 270 ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2) 275 ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1) 280 TEMP5 = ABSH*SQRT(ERK) ERR = TEMP5*(G(K)-G(KP1)) ERK = TEMP5*SIG(KP1)*GSTR(K) KNEW = K C C TEST IF ORDER SHOULD BE LOWERED C IF(KM2)299,290,285 285 IF(MAX(ERKM1,ERKM2) .LE. ERK) KNEW = KM1 GO TO 299 290 IF(ERKM1 .LE. 0.5*ERK) KNEW = KM1 C C TEST IF STEP SUCCESSFUL C 299 IF(ERR .LE. EPS) GO TO 400 C *** END BLOCK 2 *** C C *** BEGIN BLOCK 3 *** C THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), PSI(*) . C IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE. IF STEP FAILS MORE C THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE ERROR C TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE C PRECISION. C *** C C RESTORE X, PHI(*,*) AND PSI(*) C PHASE1 = .FALSE. X = XOLD DO 310 I = 1,K TEMP1 = 1.0/BETA(I) IP1 = I+1 DO 305 L = 1,NEQN 305 PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1)) 310 CONTINUE IF(K .LT. 2) GO TO 320 DO 315 I = 2,K 315 PSI(I-1) = PSI(I) - H C C ON THIRD FAILURE, SET ORDER TO ONE. THEREAFTER, USE OPTIMAL STEP C SIZE C 320 IFAIL = IFAIL + 1 TEMP2 = 0.5 IF(IFAIL - 3) 335,330,325 325 IF(P5EPS .LT. 0.25*ERK) TEMP2 = SQRT(P5EPS/ERK) 330 KNEW = 1 335 H = TEMP2*H K = KNEW NS = 0 IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340 CRASH = .TRUE. H = SIGN(FOURU*ABS(X),H) EPS = EPS + EPS RETURN 340 GO TO 100 C *** END BLOCK 3 *** C C *** BEGIN BLOCK 4 *** C THE STEP IS SUCCESSFUL. CORRECT THE PREDICTED SOLUTION, EVALUATE C THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE C DIFFERENCES. DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP. C *** 400 KOLD = K HOLD = H C C CORRECT AND EVALUATE C TEMP1 = H*G(KP1) IF(NORND) GO TO 410 DO 405 L = 1,NEQN TEMP3 = Y(L) RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16) Y(L) = P(L) + RHO PHI(L,15) = (Y(L) - P(L)) - RHO 405 P(L) = TEMP3 GO TO 420 410 DO 415 L = 1,NEQN TEMP3 = Y(L) Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1)) 415 P(L) = TEMP3 420 CALL F(X,Y,YP,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,IFPWA1, $ IFPC1,NEQN-1,IFPC2,PAR,IPAR) IF (IFPC2 .GT. 0) RETURN C C UPDATE DIFFERENCES FOR NEXT STEP C DO 425 L = 1,NEQN PHI(L,KP1) = YP(L) - PHI(L,1) 425 PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2) DO 435 I = 1,K DO 430 L = 1,NEQN 430 PHI(L,I) = PHI(L,I) + PHI(L,KP1) 435 CONTINUE C C ESTIMATE ERROR AT ORDER K+1 UNLESS: C IN FIRST PHASE WHEN ALWAYS RAISE ORDER, C ALREADY DECIDED TO LOWER ORDER, C STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE C ERKP1 = 0.0 IF(KNEW .EQ. KM1 .OR. K .EQ. 12) PHASE1 = .FALSE. IF(PHASE1) GO TO 450 IF(KNEW .EQ. KM1) GO TO 455 IF(KP1 .GT. NS) GO TO 460 DO 440 L = 1,NEQN 440 ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2 ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1) C C USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER C FOR NEXT STEP C IF(K .GT. 1) GO TO 445 IF(ERKP1 .GE. 0.5*ERK) GO TO 460 GO TO 450 445 IF(ERKM1 .LE. MIN(ERK,ERKP1)) GO TO 455 IF(ERKP1 .GE. ERK .OR. K .EQ. 12) GO TO 460 C C HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE C BEEN LOWERED IN BLOCK 2. THUS ORDER IS TO BE RAISED C C RAISE ORDER C 450 K = KP1 ERK = ERKP1 GO TO 460 C C LOWER ORDER C 455 K = KM1 ERK = ERKM1 C C WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP C 460 HNEW = H + H IF(PHASE1) GO TO 465 IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465 HNEW = H IF(P5EPS .GE. ERK) GO TO 465 TEMP2 = K+1 R = (P5EPS/ERK)**(1.0/TEMP2) HNEW = ABSH*MAX(0.5D0,MIN(0.9D0,R)) HNEW = SIGN(MAX(HNEW,FOURU*ABS(X)),H) 465 H = HNEW RETURN C *** END BLOCK 4 *** END