GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
stode.f
Go to the documentation of this file.
1  SUBROUTINE stode (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
2  1 wm, iwm, f, jac, pjac, slvs, ierr)
3 CLLL. OPTIMIZE
4  EXTERNAL f, jac, pjac, slvs
5  INTEGER neq, nyh, iwm
6  INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
7  1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
8  2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9  INTEGER i, i1, iredo, iret, j, jb, m, ncf, newq
10  DOUBLE PRECISION y, yh, yh1, ewt, savf, acor, wm
11  DOUBLE PRECISION conit, crate, el, elco, hold, rmax, tesco,
12  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
13  DOUBLE PRECISION dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
14  1 r, rh, rhdn, rhsm, rhup, told, vnorm
15  dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
16  1 acor(*), wm(*), iwm(*)
17  COMMON /ls0001/ conit, crate, el(13), elco(13,12),
18  1 hold, rmax, tesco(3,12),
19  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14),
20  3 ialth, ipup, lmax, meo, nqnyh, nslp,
21  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
22  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
23 C-----------------------------------------------------------------------
24 C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
25 C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
26 C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD
27 C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT
28 C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE.
29 C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES..
30 C
31 C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND
32 C PASSED AS THE NEQ ARGUMENT IN ALL CALLS TO F AND JAC.
33 C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN
34 C ALL CALLS TO F AND JAC.
35 C YH = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
36 C AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE
37 C LMAX = MAXORD + 1. YH(I,J+1) CONTAINS THE APPROXIMATE
38 C J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J)
39 C (J = 0,1,...,NQ). ON ENTRY FOR THE FIRST STEP, THE FIRST
40 C TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES.
41 C NYH = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH.
42 C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH.
43 C EWT = AN ARRAY OF LENGTH N CONTAINING MULTIPLICATIVE WEIGHTS
44 C FOR LOCAL ERROR MEASUREMENTS. LOCAL ERRORS IN Y(I) ARE
45 C COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS.
46 C SAVF = AN ARRAY OF WORKING STORAGE, OF LENGTH N.
47 C ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1
48 C AND MAXORD .LT. THE CURRENT ORDER NQ.
49 C ACOR = A WORK ARRAY OF LENGTH N, USED FOR THE ACCUMULATED
50 C CORRECTIONS. ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS
51 C THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I).
52 C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX
53 C OPERATIONS IN CHORD ITERATION (MITER .NE. 0).
54 C PJAC = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX
55 C AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED.
56 C SLVS = NAME OF ROUTINE TO SOLVE LINEAR SYSTEM IN CHORD ITERATION.
57 C CCMAX = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED.
58 C H = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
59 C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
60 C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
61 C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
62 C HMIN = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED.
63 C HMXI = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED.
64 C HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX.
65 C HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT
66 C TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED.
67 C TN = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN.
68 C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING
69 C VALUES AND MEANINGS..
70 C 0 PERFORM THE FIRST STEP.
71 C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST.
72 C -1 TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD,
73 C N, METH, MITER, AND/OR MATRIX PARAMETERS.
74 C -2 TAKE THE NEXT STEP WITH A NEW VALUE OF H,
75 C BUT WITH OTHER INPUTS UNCHANGED.
76 C ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION.
77 C KFLAG = A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
78 C 0 THE STEP WAS SUCCESFUL.
79 C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED.
80 C -2 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED.
81 C -3 FATAL ERROR IN PJAC OR SLVS.
82 C A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER
83 C ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED.
84 C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND
85 C THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST
86 C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
87 C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED.
88 C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED.
89 C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0).
90 C MXNCF = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED.
91 C METH/MITER = THE METHOD FLAGS. SEE DESCRIPTION IN DRIVER.
92 C N = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
93 C IERR = ERROR FLAG FROM USER-SUPPLIED FUNCTION
94 C-----------------------------------------------------------------------
95  kflag = 0
96  told = tn
97  ncf = 0
98  ierpj = 0
99  iersl = 0
100  jcur = 0
101  icf = 0
102  delp = 0.0d0
103  IF (jstart .GT. 0) go to 200
104  IF (jstart .EQ. -1) go to 100
105  IF (jstart .EQ. -2) go to 160
106 C-----------------------------------------------------------------------
107 C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
108 C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
109 C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
110 C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE
111 C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
112 C FOR THE NEXT INCREASE.
113 C-----------------------------------------------------------------------
114  lmax = maxord + 1
115  nq = 1
116  l = 2
117  ialth = 2
118  rmax = 10000.0d0
119  rc = 0.0d0
120  el0 = 1.0d0
121  crate = 0.7d0
122  hold = h
123  meo = meth
124  nslp = 0
125  ipup = miter
126  iret = 3
127  go to 140
128 C-----------------------------------------------------------------------
129 C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
130 C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
131 C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
132 C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
133 C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET
134 C THE COEFFICIENTS OF THE METHOD.
135 C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
136 C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
137 C IF H IS TO BE CHANGED, YH MUST BE RESCALED.
138 C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
139 C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
140 C-----------------------------------------------------------------------
141  100 ipup = miter
142  lmax = maxord + 1
143  IF (ialth .EQ. 1) ialth = 2
144  IF (meth .EQ. meo) go to 110
145  CALL cfode(meth, elco, tesco)
146  meo = meth
147  IF (nq .GT. maxord) go to 120
148  ialth = l
149  iret = 1
150  go to 150
151  110 IF (nq .LE. maxord) go to 160
152  120 nq = maxord
153  l = lmax
154  DO 125 i = 1,l
155  125 el(i) = elco(i,nq)
156  nqnyh = nq*nyh
157  rc = rc*el(1)/el0
158  el0 = el(1)
159  conit = 0.5d0/dble(nq+2)
160  ddn = vnorm(n, savf, ewt)/tesco(1,l)
161  exdn = 1.0d0/dble(l)
162  rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
163  rh = dmin1(rhdn,1.0d0)
164  iredo = 3
165  IF (h .EQ. hold) go to 170
166  rh = dmin1(rh,dabs(h/hold))
167  h = hold
168  go to 175
169 C-----------------------------------------------------------------------
170 C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
171 C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
172 C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
173 C-----------------------------------------------------------------------
174  140 CALL cfode(meth, elco, tesco)
175  150 DO 155 i = 1,l
176  155 el(i) = elco(i,nq)
177  nqnyh = nq*nyh
178  rc = rc*el(1)/el0
179  el0 = el(1)
180  conit = 0.5d0/dble(nq+2)
181  go to(160, 170, 200), iret
182 C-----------------------------------------------------------------------
183 C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
184 C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO
185 C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
186 C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
187 C-----------------------------------------------------------------------
188  160 IF (h .EQ. hold) go to 200
189  rh = h/hold
190  h = hold
191  iredo = 3
192  go to 175
193  170 rh = dmax1(rh,hmin/dabs(h))
194  175 rh = dmin1(rh,rmax)
195  rh = rh/dmax1(1.0d0,dabs(h)*hmxi*rh)
196  r = 1.0d0
197  DO 180 j = 2,l
198  r = r*rh
199  DO 180 i = 1,n
200  180 yh(i,j) = yh(i,j)*r
201  h = h*rh
202  rc = rc*rh
203  ialth = l
204  IF (iredo .EQ. 0) go to 690
205 C-----------------------------------------------------------------------
206 C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
207 C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
208 C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1).
209 C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER
210 C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
211 C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS.
212 C-----------------------------------------------------------------------
213  200 IF (dabs(rc-1.0d0) .GT. ccmax) ipup = miter
214  IF (nst .GE. nslp+msbp) ipup = miter
215  tn = tn + h
216  i1 = nqnyh + 1
217  DO 215 jb = 1,nq
218  i1 = i1 - nyh
219 CDIR$ IVDEP
220  DO 210 i = i1,nqnyh
221  210 yh1(i) = yh1(i) + yh1(i+nyh)
222  215 CONTINUE
223 C-----------------------------------------------------------------------
224 C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS
225 C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR
226 C WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE
227 C VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP.
228 C-----------------------------------------------------------------------
229  220 m = 0
230  DO 230 i = 1,n
231  230 y(i) = yh(i,1)
232  ierr = 0
233  CALL f(neq, tn, y, savf, ierr)
234  IF (ierr .LT. 0) RETURN
235  nfe = nfe + 1
236  IF (ipup .LE. 0) go to 250
237 C-----------------------------------------------------------------------
238 C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND
239 C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET
240 C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
241 C-----------------------------------------------------------------------
242  ierr = 0
243  CALL pjac(neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac,
244  1 ierr)
245  IF (ierr .LT. 0) RETURN
246  ipup = 0
247  rc = 1.0d0
248  nslp = nst
249  crate = 0.7d0
250  IF (ierpj .NE. 0) go to 430
251  250 DO 260 i = 1,n
252  260 acor(i) = 0.0d0
253  270 IF (miter .NE. 0) go to 350
254 C-----------------------------------------------------------------------
255 C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
256 C THE RESULT OF THE LAST FUNCTION EVALUATION.
257 C-----------------------------------------------------------------------
258  DO 290 i = 1,n
259  savf(i) = h*savf(i) - yh(i,2)
260  290 y(i) = savf(i) - acor(i)
261  del = vnorm(n, y, ewt)
262  DO 300 i = 1,n
263  y(i) = yh(i,1) + el(1)*savf(i)
264  300 acor(i) = savf(i)
265  go to 400
266 C-----------------------------------------------------------------------
267 C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
268 C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
269 C P AS COEFFICIENT MATRIX.
270 C-----------------------------------------------------------------------
271  350 DO 360 i = 1,n
272  360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
273  CALL slvs(wm, iwm, y, savf)
274  IF (iersl .LT. 0) go to 430
275  IF (iersl .GT. 0) go to 410
276  del = vnorm(n, y, ewt)
277  DO 380 i = 1,n
278  acor(i) = acor(i) + y(i)
279  380 y(i) = yh(i,1) + el(1)*acor(i)
280 C-----------------------------------------------------------------------
281 C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
282 C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
283 C-----------------------------------------------------------------------
284  400 IF (m .NE. 0) crate = dmax1(0.2d0*crate,del/delp)
285  dcon = del*dmin1(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit)
286  IF (dcon .LE. 1.0d0) go to 450
287  m = m + 1
288  IF (m .EQ. maxcor) go to 410
289  IF (m .GE. 2 .AND. del .GT. 2.0d0*delp) go to 410
290  delp = del
291  ierr = 0
292  CALL f(neq, tn, y, savf, ierr)
293  IF (ierr .LT. 0) RETURN
294  nfe = nfe + 1
295  go to 270
296 C-----------------------------------------------------------------------
297 C THE CORRECTOR ITERATION FAILED TO CONVERGE.
298 C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR
299 C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
300 C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE
301 C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
302 C-----------------------------------------------------------------------
303  410 IF (miter .EQ. 0 .OR. jcur .EQ. 1) go to 430
304  icf = 1
305  ipup = miter
306  go to 220
307  430 icf = 2
308  ncf = ncf + 1
309  rmax = 2.0d0
310  tn = told
311  i1 = nqnyh + 1
312  DO 445 jb = 1,nq
313  i1 = i1 - nyh
314 CDIR$ IVDEP
315  DO 440 i = i1,nqnyh
316  440 yh1(i) = yh1(i) - yh1(i+nyh)
317  445 CONTINUE
318  IF (ierpj .LT. 0 .OR. iersl .LT. 0) go to 680
319  IF (dabs(h) .LE. hmin*1.00001d0) go to 670
320  IF (ncf .EQ. mxncf) go to 670
321  rh = 0.25d0
322  ipup = miter
323  iredo = 1
324  go to 170
325 C-----------------------------------------------------------------------
326 C THE CORRECTOR HAS CONVERGED. JCUR IS SET TO 0
327 C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER.
328 C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
329 C IF IT FAILS.
330 C-----------------------------------------------------------------------
331  450 jcur = 0
332  IF (m .EQ. 0) dsm = del/tesco(2,nq)
333  IF (m .GT. 0) dsm = vnorm(n, acor, ewt)/tesco(2,nq)
334  IF (dsm .GT. 1.0d0) go to 500
335 C-----------------------------------------------------------------------
336 C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
337 C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1.
338 C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
339 C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
340 C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
341 C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A
342 C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT
343 C TESTING FOR THAT MANY STEPS.
344 C-----------------------------------------------------------------------
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 700
356  IF (l .EQ. lmax) go to 700
357  DO 490 i = 1,n
358  490 yh(i,lmax) = acor(i)
359  go to 700
360 C-----------------------------------------------------------------------
361 C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
362 C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
363 C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
364 C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE
365 C BY A FACTOR OF 0.2 OR LESS.
366 C-----------------------------------------------------------------------
367  500 kflag = kflag - 1
368  tn = told
369  i1 = nqnyh + 1
370  DO 515 jb = 1,nq
371  i1 = i1 - nyh
372 CDIR$ IVDEP
373  DO 510 i = i1,nqnyh
374  510 yh1(i) = yh1(i) - yh1(i+nyh)
375  515 CONTINUE
376  rmax = 2.0d0
377  IF (dabs(h) .LE. hmin*1.00001d0) go to 660
378  IF (kflag .LE. -3) go to 640
379  iredo = 2
380  rhup = 0.0d0
381  go to 540
382 C-----------------------------------------------------------------------
383 C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
384 C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
385 C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
386 C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
387 C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
388 C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
389 C ADDITIONAL SCALED DERIVATIVE.
390 C-----------------------------------------------------------------------
391  520 rhup = 0.0d0
392  IF (l .EQ. lmax) go to 540
393  DO 530 i = 1,n
394  530 savf(i) = acor(i) - yh(i,lmax)
395  dup = vnorm(n, savf, ewt)/tesco(3,nq)
396  exup = 1.0d0/dble(l+1)
397  rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
398  540 exsm = 1.0d0/dble(l)
399  rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
400  rhdn = 0.0d0
401  IF (nq .EQ. 1) go to 560
402  ddn = vnorm(n, yh(1,l), ewt)/tesco(1,nq)
403  exdn = 1.0d0/dble(nq)
404  rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
405  560 IF (rhsm .GE. rhup) go to 570
406  IF (rhup .GT. rhdn) go to 590
407  go to 580
408  570 IF (rhsm .LT. rhdn) go to 580
409  newq = nq
410  rh = rhsm
411  go to 620
412  580 newq = nq - 1
413  rh = rhdn
414  IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
415  go to 620
416  590 newq = l
417  rh = rhup
418  IF (rh .LT. 1.1d0) go to 610
419  r = el(l)/dble(l)
420  DO 600 i = 1,n
421  600 yh(i,newq+1) = acor(i)*r
422  go to 630
423  610 ialth = 3
424  go to 700
425  620 IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1d0)) go to 610
426  IF (kflag .LE. -2) rh = dmin1(rh,0.2d0)
427 C-----------------------------------------------------------------------
428 C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
429 C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
430 C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
431 C-----------------------------------------------------------------------
432  IF (newq .EQ. nq) go to 170
433  630 nq = newq
434  l = nq + 1
435  iret = 2
436  go to 150
437 C-----------------------------------------------------------------------
438 C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
439 C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
440 C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
441 C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
442 C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN
443 C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
444 C UNTIL IT SUCCEEDS OR H REACHES HMIN.
445 C-----------------------------------------------------------------------
446  640 IF (kflag .EQ. -10) go to 660
447  rh = 0.1d0
448  rh = dmax1(hmin/dabs(h),rh)
449  h = h*rh
450  DO 645 i = 1,n
451  645 y(i) = yh(i,1)
452  ierr = 0
453  CALL f(neq, tn, y, savf, ierr)
454  IF (ierr .LT. 0) RETURN
455  nfe = nfe + 1
456  DO 650 i = 1,n
457  650 yh(i,2) = h*savf(i)
458  ipup = miter
459  ialth = 5
460  IF (nq .EQ. 1) go to 200
461  nq = 1
462  l = 2
463  iret = 3
464  go to 150
465 C-----------------------------------------------------------------------
466 C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD
467 C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
468 C-----------------------------------------------------------------------
469  660 kflag = -1
470  go to 720
471  670 kflag = -2
472  go to 720
473  680 kflag = -3
474  go to 720
475  690 rmax = 10.0d0
476  700 r = 1.0d0/tesco(2,nqu)
477  DO 710 i = 1,n
478  710 acor(i) = acor(i)*r
479  720 hold = h
480  jstart = 1
481  RETURN
482 C----------------------- END OF SUBROUTINE STODE -----------------------
483  END