1*DECK STOD
2      SUBROUTINE STOD (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM,
3     +   F, JAC, RPAR, IPAR)
4C***BEGIN PROLOGUE  STOD
5C***SUBSIDIARY
6C***PURPOSE  Subsidiary to DEBDF
7C***LIBRARY   SLATEC
8C***TYPE      SINGLE PRECISION (STOD-S, DSTOD-D)
9C***AUTHOR  Watts, H. A., (SNLA)
10C***DESCRIPTION
11C
12C   STOD integrates a system of first order odes over one step in the
13C   integrator package DEBDF.
14C ----------------------------------------------------------------------
15C STOD  performs one step of the integration of an initial value
16C problem for a system of ordinary differential equations.
17C Note.. STOD  is independent of the value of the iteration method
18C indicator MITER, when this is .NE. 0, and hence is independent
19C of the type of chord method used, or the Jacobian structure.
20C Communication with STOD  is done with the following variables..
21C
22C Y      = An array of length .GE. n used as the Y argument in
23C          all calls to F and JAC.
24C NEQ    = Integer array containing problem size in NEQ(1), and
25C          passed as the NEQ argument in all calls to F and JAC.
26C YH     = An NYH by LMAX array containing the dependent variables
27C          and their approximate scaled derivatives, where
28C          LMAX = MAXORD + 1.  YH(I,J+1) contains the approximate
29C          J-th derivative of Y(I), scaled by H**J/Factorial(j)
30C          (J = 0,1,...,NQ).  On entry for the first step, the first
31C          two columns of YH must be set from the initial values.
32C NYH    = A constant integer .GE. N, the first dimension of YH.
33C YH1    = A one-dimensional array occupying the same space as YH.
34C EWT    = An array of N elements with which the estimated local
35C          errors in YH are compared.
36C SAVF   = An array of working storage, of length N.
37C ACOR   = A work array of length N, used for the accumulated
38C          corrections.  On a successful return, ACOR(I) contains
39C          the estimated one-step local error in Y(I).
40C WM,IWM = Real and integer work arrays associated with matrix
41C          operations in chord iteration (MITER .NE. 0).
42C PJAC   = Name of routine to evaluate and preprocess Jacobian matrix
43C          if a chord method is being used.
44C SLVS   = Name of routine to solve linear system in chord iteration.
45C H      = The step size to be attempted on the next step.
46C          H is altered by the error control algorithm during the
47C          problem.  H can be either positive or negative, but its
48C          sign must remain constant throughout the problem.
49C HMIN   = The minimum absolute value of the step size H to be used.
50C HMXI   = Inverse of the maximum absolute value of H to be used.
51C          HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
52C          HMIN and HMXI may be changed at any time, but will not
53C          take effect until the next change of H is considered.
54C TN     = The independent variable. TN is updated on each step taken.
55C JSTART = An integer used for input only, with the following
56C          values and meanings..
57C               0  Perform the first step.
58C           .GT.0  Take a new step continuing from the last.
59C              -1  Take the next step with a new value of H, MAXORD,
60C                    N, METH, MITER, and/or matrix parameters.
61C              -2  Take the next step with a new value of H,
62C                    but with other inputs unchanged.
63C          On return, JSTART is set to 1 to facilitate continuation.
64C KFLAG  = a completion code with the following meanings..
65C               0  The step was successful.
66C              -1  The requested error could not be achieved.
67C              -2  Corrector convergence could not be achieved.
68C          A return with KFLAG = -1 or -2 means either
69C          ABS(H) = HMIN or 10 consecutive failures occurred.
70C          On a return with KFLAG negative, the values of TN and
71C          the YH array are as of the beginning of the last
72C          step, and H is the last step size attempted.
73C MAXORD = The maximum order of integration method to be allowed.
74C METH/MITER = The method flags.  See description in driver.
75C N      = The number of first-order differential equations.
76C ----------------------------------------------------------------------
77C
78C***SEE ALSO  DEBDF
79C***ROUTINES CALLED  CFOD, PJAC, SLVS, VNWRMS
80C***COMMON BLOCKS    DEBDF1
81C***REVISION HISTORY  (YYMMDD)
82C   800901  DATE WRITTEN
83C   890531  Changed all specific intrinsics to generic.  (WRB)
84C   891214  Prologue converted to Version 4.0 format.  (BAB)
85C   900328  Added TYPE section.  (WRB)
86C   910722  Updated AUTHOR section.  (ALS)
87C   920422  Changed DIMENSION statement.  (WRB)
88C***END PROLOGUE  STOD
89      EXTERNAL F, JAC
90C
91CLLL. OPTIMIZE
92      INTEGER NEQ, NYH, IWM, I, I1, IALTH, IER, IOWND, IREDO, IRET,
93     1   IPUP, J, JB, JSTART, KFLAG, L, LMAX, M, MAXORD, MEO, METH,
94     2   MITER, N, NCF, NEWQ, NFE, NJE, NQ, NQNYH, NQU, NST, NSTEPJ
95      REAL Y, YH, YH1, EWT, SAVF, ACOR, WM,
96     1   ROWND, CONIT, CRATE, EL, ELCO, HOLD, RC, RMAX, TESCO,
97     2   EL0, H, HMIN, HMXI, HU, TN, UROUND,
98     3   DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
99     4   R, RH, RHDN, RHSM, RHUP, TOLD, VNWRMS
100      DIMENSION         Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
101     1   ACOR(*), WM(*), IWM(*), RPAR(*), IPAR(*)
102      COMMON /DEBDF1/ ROWND, CONIT, CRATE, EL(13), ELCO(13,12),
103     1   HOLD, RC, RMAX, TESCO(3,12),
104     2   EL0, H, HMIN, HMXI, HU, TN, UROUND, IOWND(7), KSTEPS, IOD(6),
105     3   IALTH, IPUP, LMAX, MEO, NQNYH, NSTEPJ,
106     4   IER, JSTART, KFLAG, L, METH, MITER, MAXORD, N, NQ, NST, NFE,
107     5   NJE, NQU
108C
109C
110C***FIRST EXECUTABLE STATEMENT  STOD
111      KFLAG = 0
112      TOLD = TN
113      NCF = 0
114      IF (JSTART .GT. 0) GO TO 200
115      IF (JSTART .EQ. -1) GO TO 100
116      IF (JSTART .EQ. -2) GO TO 160
117C-----------------------------------------------------------------------
118C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
119C INITIALIZED.  RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
120C IN A SINGLE STEP.  IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
121C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10.  IF A FAILURE
122C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
123C FOR THE NEXT INCREASE.
124C-----------------------------------------------------------------------
125      LMAX = MAXORD + 1
126      NQ = 1
127      L = 2
128      IALTH = 2
129      RMAX = 10000.0E0
130      RC = 0.0E0
131      EL0 = 1.0E0
132      CRATE = 0.7E0
133      DELP = 0.0E0
134      HOLD = H
135      MEO = METH
136      NSTEPJ = 0
137      IRET = 3
138      GO TO 140
139C-----------------------------------------------------------------------
140C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
141C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
142C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
143C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
144C IF THE CALLER HAS CHANGED METH, CFOD  IS CALLED TO RESET
145C THE COEFFICIENTS OF THE METHOD.
146C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
147C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
148C IF H IS TO BE CHANGED, YH MUST BE RESCALED.
149C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
150C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
151C-----------------------------------------------------------------------
152 100  IPUP = MITER
153      LMAX = MAXORD + 1
154      IF (IALTH .EQ. 1) IALTH = 2
155      IF (METH .EQ. MEO) GO TO 110
156      CALL CFOD  (METH, ELCO, TESCO)
157      MEO = METH
158      IF (NQ .GT. MAXORD) GO TO 120
159      IALTH = L
160      IRET = 1
161      GO TO 150
162 110  IF (NQ .LE. MAXORD) GO TO 160
163 120  NQ = MAXORD
164      L = LMAX
165      DO 125 I = 1,L
166 125    EL(I) = ELCO(I,NQ)
167      NQNYH = NQ*NYH
168      RC = RC*EL(1)/EL0
169      EL0 = EL(1)
170      CONIT = 0.5E0/(NQ+2)
171      DDN = VNWRMS (N, SAVF, EWT)/TESCO(1,L)
172      EXDN = 1.0E0/L
173      RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0)
174      RH = MIN(RHDN,1.0E0)
175      IREDO = 3
176      IF (H .EQ. HOLD) GO TO 170
177      RH = MIN(RH,ABS(H/HOLD))
178      H = HOLD
179      GO TO 175
180C-----------------------------------------------------------------------
181C CFOD  IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
182C CURRENT METH.  THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
183C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
184C-----------------------------------------------------------------------
185 140  CALL CFOD  (METH, ELCO, TESCO)
186 150  DO 155 I = 1,L
187 155    EL(I) = ELCO(I,NQ)
188      NQNYH = NQ*NYH
189      RC = RC*EL(1)/EL0
190      EL0 = EL(1)
191      CONIT = 0.5E0/(NQ+2)
192      GO TO (160, 170, 200), IRET
193C-----------------------------------------------------------------------
194C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
195C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED.  IALTH IS SET TO
196C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
197C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
198C-----------------------------------------------------------------------
199 160  IF (H .EQ. HOLD) GO TO 200
200      RH = H/HOLD
201      H = HOLD
202      IREDO = 3
203      GO TO 175
204 170  RH = MAX(RH,HMIN/ABS(H))
205 175  RH = MIN(RH,RMAX)
206      RH = RH/MAX(1.0E0,ABS(H)*HMXI*RH)
207      R = 1.0E0
208      DO 180 J = 2,L
209        R = R*RH
210        DO 180 I = 1,N
211 180      YH(I,J) = YH(I,J)*R
212      H = H*RH
213      RC = RC*RH
214      IALTH = L
215      IF (IREDO .EQ. 0) GO TO 680
216C-----------------------------------------------------------------------
217C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
218C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
219C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT  H*EL(1).
220C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER
221C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
222C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY 20-TH STEP.
223C-----------------------------------------------------------------------
224 200  IF (ABS(RC-1.0E0) .GT. 0.3E0) IPUP = MITER
225      IF (NST .GE. NSTEPJ+20) IPUP = MITER
226      TN = TN + H
227      I1 = NQNYH + 1
228      DO 215 JB = 1,NQ
229        I1 = I1 - NYH
230        DO 210 I = I1,NQNYH
231 210      YH1(I) = YH1(I) + YH1(I+NYH)
232 215    CONTINUE
233      KSTEPS = KSTEPS + 1
234C-----------------------------------------------------------------------
235C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN.  A CONVERGENCE TEST IS
236C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR
237C WEIGHT VECTOR EWT.  THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE
238C VECTOR ACOR(I).  THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP.
239C-----------------------------------------------------------------------
240 220  M = 0
241      DO 230 I = 1,N
242 230    Y(I) = YH(I,1)
243      CALL F (TN, Y, SAVF, RPAR, IPAR)
244      NFE = NFE + 1
245      IF (IPUP .LE. 0) GO TO 250
246C-----------------------------------------------------------------------
247C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND
248C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION.  IPUP IS SET
249C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
250C-----------------------------------------------------------------------
251      IPUP = 0
252      RC = 1.0E0
253      NSTEPJ = NST
254      CRATE = 0.7E0
255      CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC,
256     1          RPAR, IPAR)
257      IF (IER .NE. 0) GO TO 430
258 250  DO 260 I = 1,N
259 260    ACOR(I) = 0.0E0
260 270  IF (MITER .NE. 0) GO TO 350
261C-----------------------------------------------------------------------
262C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
263C THE RESULT OF THE LAST FUNCTION EVALUATION.
264C-----------------------------------------------------------------------
265      DO 290 I = 1,N
266        SAVF(I) = H*SAVF(I) - YH(I,2)
267 290    Y(I) = SAVF(I) - ACOR(I)
268      DEL = VNWRMS (N, Y, EWT)
269      DO 300 I = 1,N
270        Y(I) = YH(I,1) + EL(1)*SAVF(I)
271 300    ACOR(I) = SAVF(I)
272      GO TO 400
273C-----------------------------------------------------------------------
274C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
275C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
276C P AS COEFFICIENT MATRIX.
277C-----------------------------------------------------------------------
278 350  DO 360 I = 1,N
279 360    Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
280      CALL SLVS (WM, IWM, Y, SAVF)
281      IF (IER .NE. 0) GO TO 410
282      DEL = VNWRMS (N, Y, EWT)
283      DO 380 I = 1,N
284        ACOR(I) = ACOR(I) + Y(I)
285 380    Y(I) = YH(I,1) + EL(1)*ACOR(I)
286C-----------------------------------------------------------------------
287C TEST FOR CONVERGENCE.  IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
288C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
289C-----------------------------------------------------------------------
290 400  IF (M .NE. 0) CRATE = MAX(0.2E0*CRATE,DEL/DELP)
291      DCON = DEL*MIN(1.0E0,1.5E0*CRATE)/(TESCO(2,NQ)*CONIT)
292      IF (DCON .LE. 1.0E0) GO TO 450
293      M = M + 1
294      IF (M .EQ. 3) GO TO 410
295      IF (M .GE. 2 .AND. DEL .GT. 2.0E0*DELP) GO TO 410
296      DELP = DEL
297      CALL F (TN, Y, SAVF, RPAR, IPAR)
298      NFE = NFE + 1
299      GO TO 270
300C-----------------------------------------------------------------------
301C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES.
302C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR
303C THE NEXT TRY.  OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
304C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE.  IF H CANNOT BE
305C REDUCED OR 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
306C-----------------------------------------------------------------------
307 410  IF (IPUP .EQ. 0) GO TO 430
308      IPUP = MITER
309      GO TO 220
310 430  TN = TOLD
311      NCF = NCF + 1
312      RMAX = 2.0E0
313      I1 = NQNYH + 1
314      DO 445 JB = 1,NQ
315        I1 = I1 - NYH
316        DO 440 I = I1,NQNYH
317 440      YH1(I) = YH1(I) - YH1(I+NYH)
318 445    CONTINUE
319      IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 670
320      IF (NCF .EQ. 10) GO TO 670
321      RH = 0.25E0
322      IPUP = MITER
323      IREDO = 1
324      GO TO 170
325C-----------------------------------------------------------------------
326C THE CORRECTOR HAS CONVERGED.  IPUP IS SET TO -1 IF MITER .NE. 0,
327C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER.
328C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
329C IF IT FAILS.
330C-----------------------------------------------------------------------
331 450  IF (MITER .NE. 0) IPUP = -1
332      IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
333      IF (M .GT. 0) DSM = VNWRMS (N, ACOR, EWT)/TESCO(2,NQ)
334      IF (DSM .GT. 1.0E0) GO TO 500
335C-----------------------------------------------------------------------
336C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
337C CONSIDER CHANGING H IF IALTH = 1.  OTHERWISE DECREASE IALTH BY 1.
338C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
339C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
340C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
341C BY ONE IS CONSIDERED ALSO.  A CHANGE IN H IS MADE ONLY IF IT IS BY A
342C FACTOR OF AT LEAST 1.1.  IF NOT, IALTH IS SET TO 3 TO PREVENT
343C TESTING FOR THAT MANY STEPS.
344C-----------------------------------------------------------------------
345      KFLAG = 0
346      IREDO = 0
347      NST = NST + 1
348      HU = H
349      NQU = NQ
350      DO 470 J = 1,L
351        DO 470 I = 1,N
352 470      YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
353      IALTH = IALTH - 1
354      IF (IALTH .EQ. 0) GO TO 520
355      IF (IALTH .GT. 1) GO TO 690
356      IF (L .EQ. LMAX) GO TO 690
357      DO 490 I = 1,N
358 490    YH(I,LMAX) = ACOR(I)
359      GO TO 690
360C-----------------------------------------------------------------------
361C THE ERROR TEST FAILED.  KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
362C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
363C TO TRY THE STEP AGAIN.  COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
364C ONE LOWER ORDER.  AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE
365C BY A FACTOR OF 0.2 OR LESS.
366C-----------------------------------------------------------------------
367 500  KFLAG = KFLAG - 1
368      TN = TOLD
369      I1 = NQNYH + 1
370      DO 515 JB = 1,NQ
371        I1 = I1 - NYH
372        DO 510 I = I1,NQNYH
373 510      YH1(I) = YH1(I) - YH1(I+NYH)
374 515    CONTINUE
375      RMAX = 2.0E0
376      IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 660
377      IF (KFLAG .LE. -3) GO TO 640
378      IREDO = 2
379      RHUP = 0.0E0
380      GO TO 540
381C-----------------------------------------------------------------------
382C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
383C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
384C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
385C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
386C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
387C ACCORDINGLY.  IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
388C ADDITIONAL SCALED DERIVATIVE.
389C-----------------------------------------------------------------------
390 520  RHUP = 0.0E0
391      IF (L .EQ. LMAX) GO TO 540
392      DO 530 I = 1,N
393 530    SAVF(I) = ACOR(I) - YH(I,LMAX)
394      DUP = VNWRMS (N, SAVF, EWT)/TESCO(3,NQ)
395      EXUP = 1.0E0/(L+1)
396      RHUP = 1.0E0/(1.4E0*DUP**EXUP + 0.0000014E0)
397 540  EXSM = 1.0E0/L
398      RHSM = 1.0E0/(1.2E0*DSM**EXSM + 0.0000012E0)
399      RHDN = 0.0E0
400      IF (NQ .EQ. 1) GO TO 560
401      DDN = VNWRMS (N, YH(1,L), EWT)/TESCO(1,NQ)
402      EXDN = 1.0E0/NQ
403      RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0)
404 560  IF (RHSM .GE. RHUP) GO TO 570
405      IF (RHUP .GT. RHDN) GO TO 590
406      GO TO 580
407 570  IF (RHSM .LT. RHDN) GO TO 580
408      NEWQ = NQ
409      RH = RHSM
410      GO TO 620
411 580  NEWQ = NQ - 1
412      RH = RHDN
413      IF (KFLAG .LT. 0 .AND. RH .GT. 1.0E0) RH = 1.0E0
414      GO TO 620
415 590  NEWQ = L
416      RH = RHUP
417      IF (RH .LT. 1.1E0) GO TO 610
418      R = EL(L)/L
419      DO 600 I = 1,N
420 600    YH(I,NEWQ+1) = ACOR(I)*R
421      GO TO 630
422 610  IALTH = 3
423      GO TO 690
424 620  IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1E0)) GO TO 610
425      IF (KFLAG .LE. -2) RH = MIN(RH,0.2E0)
426C-----------------------------------------------------------------------
427C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
428C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
429C THEN EXIT FROM 680 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
430C-----------------------------------------------------------------------
431      IF (NEWQ .EQ. NQ) GO TO 170
432 630  NQ = NEWQ
433      L = NQ + 1
434      IRET = 2
435      GO TO 150
436C-----------------------------------------------------------------------
437C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURRED.
438C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
439C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
440C YH ARRAY HAVE ERRORS OF THE WRONG ORDER.  HENCE THE FIRST
441C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1.  THEN
442C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
443C UNTIL IT SUCCEEDS OR H REACHES HMIN.
444C-----------------------------------------------------------------------
445 640  IF (KFLAG .EQ. -10) GO TO 660
446      RH = 0.1E0
447      RH = MAX(HMIN/ABS(H),RH)
448      H = H*RH
449      DO 645 I = 1,N
450 645    Y(I) = YH(I,1)
451      CALL F (TN, Y, SAVF, RPAR, IPAR)
452      NFE = NFE + 1
453      DO 650 I = 1,N
454 650    YH(I,2) = H*SAVF(I)
455      IPUP = MITER
456      IALTH = 5
457      IF (NQ .EQ. 1) GO TO 200
458      NQ = 1
459      L = 2
460      IRET = 3
461      GO TO 150
462C-----------------------------------------------------------------------
463C ALL RETURNS ARE MADE THROUGH THIS SECTION.  H IS SAVED IN HOLD
464C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
465C-----------------------------------------------------------------------
466 660  KFLAG = -1
467      GO TO 700
468 670  KFLAG = -2
469      GO TO 700
470 680  RMAX = 10.0E0
471 690  R = 1.0E0/TESCO(2,NQU)
472      DO 695 I = 1,N
473 695    ACOR(I) = ACOR(I)*R
474 700  HOLD = H
475      JSTART = 1
476      RETURN
477C----------------------- END OF SUBROUTINE STOD  -----------------------
478      END
479