share/odepack/fortran/dstode.f

   1 *DECK DSTODE
   2       SUBROUTINE DSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
   3      1   WM, IWM, F, JAC, PJAC, SLVS)
   4 C***BEGIN PROLOGUE  DSTODE
   5 C***SUBSIDIARY
   6 C***PURPOSE  Performs one step of an ODEPACK integration.
   7 C***TYPE      DOUBLE PRECISION (SSTODE-S, DSTODE-D)
   8 C***AUTHOR  Hindmarsh, Alan C., (LLNL)
   9 C***DESCRIPTION
  10 C
  11 C  DSTODE performs one step of the integration of an initial value
  12 C  problem for a system of ordinary differential equations.
  13 C  Note:  DSTODE is independent of the value of the iteration method
  14 C  indicator MITER, when this is .ne. 0, and hence is independent
  15 C  of the type of chord method used, or the Jacobian structure.
  16 C  Communication with DSTODE is done with the following variables:
  17 C
  18 C  NEQ    = integer array containing problem size in NEQ(1), and
  19 C           passed as the NEQ argument in all calls to F and JAC.
  20 C  Y      = an array of length .ge. N used as the Y argument in
  21 C           all calls to F and JAC.
  22 C  YH     = an NYH by LMAX array containing the dependent variables
  23 C           and their approximate scaled derivatives, where
  24 C           LMAX = MAXORD + 1.  YH(i,j+1) contains the approximate
  25 C           j-th derivative of y(i), scaled by h**j/factorial(j)
  26 C           (j = 0,1,...,NQ).  on entry for the first step, the first
  27 C           two columns of YH must be set from the initial values.
  28 C  NYH    = a constant integer .ge. N, the first dimension of YH.
  29 C  YH1    = a one-dimensional array occupying the same space as YH.
  30 C  EWT    = an array of length N containing multiplicative weights
  31 C           for local error measurements.  Local errors in Y(i) are
  32 C           compared to 1.0/EWT(i) in various error tests.
  33 C  SAVF   = an array of working storage, of length N.
  34 C           Also used for input of YH(*,MAXORD+2) when JSTART = -1
  35 C           and MAXORD .lt. the current order NQ.
  36 C  ACOR   = a work array of length N, used for the accumulated
  37 C           corrections.  On a successful return, ACOR(i) contains
  38 C           the estimated one-step local error in Y(i).
  39 C  WM,IWM = real and integer work arrays associated with matrix
  40 C           operations in chord iteration (MITER .ne. 0).
  41 C  PJAC   = name of routine to evaluate and preprocess Jacobian matrix
  42 C           and P = I - h*el0*JAC, if a chord method is being used.
  43 C  SLVS   = name of routine to solve linear system in chord iteration.
  44 C  CCMAX  = maximum relative change in h*el0 before PJAC is called.
  45 C  H      = the step size to be attempted on the next step.
  46 C           H is altered by the error control algorithm during the
  47 C           problem.  H can be either positive or negative, but its
  48 C           sign must remain constant throughout the problem.
  49 C  HMIN   = the minimum absolute value of the step size h to be used.
  50 C  HMXI   = inverse of the maximum absolute value of h to be used.
  51 C           HMXI = 0.0 is allowed and corresponds to an infinite hmax.
  52 C           HMIN and HMXI may be changed at any time, but will not
  53 C           take effect until the next change of h is considered.
  54 C  TN     = the independent variable. TN is updated on each step taken.
  55 C  JSTART = an integer used for input only, with the following
  56 C           values and meanings:
  57 C                0  perform the first step.
  58 C            .gt.0  take a new step continuing from the last.
  59 C               -1  take the next step with a new value of H, MAXORD,
  60 C                     N, METH, MITER, and/or matrix parameters.
  61 C               -2  take the next step with a new value of H,
  62 C                     but with other inputs unchanged.
  63 C           On return, JSTART is set to 1 to facilitate continuation.
  64 C  KFLAG  = a completion code with the following meanings:
  65 C                0  the step was succesful.
  66 C               -1  the requested error could not be achieved.
  67 C               -2  corrector convergence could not be achieved.
  68 C               -3  fatal error in PJAC or SLVS.
  69 C           A return with KFLAG = -1 or -2 means either
  70 C           abs(H) = HMIN or 10 consecutive failures occurred.
  71 C           On a return with KFLAG negative, the values of TN and
  72 C           the YH array are as of the beginning of the last
  73 C           step, and H is the last step size attempted.
  74 C  MAXORD = the maximum order of integration method to be allowed.
  75 C  MAXCOR = the maximum number of corrector iterations allowed.
  76 C  MSBP   = maximum number of steps between PJAC calls (MITER .gt. 0).
  77 C  MXNCF  = maximum number of convergence failures allowed.
  78 C  METH/MITER = the method flags.  See description in driver.
  79 C  N      = the number of first-order differential equations.
  80 C  The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
  81 C  MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
  82 C
  83 C***SEE ALSO  DLSODE
  84 C***ROUTINES CALLED  DCFODE, DVNORM
  85 C***COMMON BLOCKS    DLS001
  86 C***REVISION HISTORY  (YYMMDD)
  87 C   791129  DATE WRITTEN
  88 C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
  89 C   890503  Minor cosmetic changes.  (FNF)
  90 C   930809  Renamed to allow single/double precision versions. (ACH)
  91 C   010418  Reduced size of Common block /DLS001/. (ACH)
  92 C   031105  Restored 'own' variables to Common block /DLS001/, to
  93 C           enable interrupt/restart feature. (ACH)
  94 C***END PROLOGUE  DSTODE
  95 C**End
  96       EXTERNAL F, JAC, PJAC, SLVS
  97       INTEGER NEQ, NYH, IWM
  98       DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
  99       DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
 100      1   ACOR(*), WM(*), IWM(*)
 101       INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
 102      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
 103      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
 104      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
 105       INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
 106       DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
 107      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
 108       DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
 109      1   R, RH, RHDN, RHSM, RHUP, TOLD, DVNORM
 110       COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12),
 111      1   HOLD, RMAX, TESCO(3,12),
 112      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
 113      3   IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
 114      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
 115      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
 116      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
 117 C
 118 C***FIRST EXECUTABLE STATEMENT  DSTODE
 119       KFLAG = 0
 120       TOLD = TN
 121       NCF = 0
 122       IERPJ = 0
 123       IERSL = 0
 124       JCUR = 0
 125       ICF = 0
 126       DELP = 0.0D0
 127       IF (JSTART .GT. 0) GO TO 200
 128       IF (JSTART .EQ. -1) GO TO 100
 129       IF (JSTART .EQ. -2) GO TO 160
 130 C-----------------------------------------------------------------------
 131 C On the first call, the order is set to 1, and other variables are
 132 C initialized.  RMAX is the maximum ratio by which H can be increased
 133 C in a single step.  It is initially 1.E4 to compensate for the small
 134 C initial H, but then is normally equal to 10.  If a failure
 135 C occurs (in corrector convergence or error test), RMAX is set to 2
 136 C for the next increase.
 137 C-----------------------------------------------------------------------
 138       LMAX = MAXORD + 1
 139       NQ = 1
 140       L = 2
 141       IALTH = 2
 142       RMAX = 10000.0D0
 143       RC = 0.0D0
 144       EL0 = 1.0D0
 145       CRATE = 0.7D0
 146       HOLD = H
 147       MEO = METH
 148       NSLP = 0
 149       IPUP = MITER
 150       IRET = 3
 151       GO TO 140
 152 C-----------------------------------------------------------------------
 153 C The following block handles preliminaries needed when JSTART = -1.
 154 C IPUP is set to MITER to force a matrix update.
 155 C If an order increase is about to be considered (IALTH = 1),
 156 C IALTH is reset to 2 to postpone consideration one more step.
 157 C If the caller has changed METH, DCFODE is called to reset
 158 C the coefficients of the method.
 159 C If the caller has changed MAXORD to a value less than the current
 160 C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
 161 C If H is to be changed, YH must be rescaled.
 162 C If H or METH is being changed, IALTH is reset to L = NQ + 1
 163 C to prevent further changes in H for that many steps.
 164 C-----------------------------------------------------------------------
 165  100  IPUP = MITER
 166       LMAX = MAXORD + 1
 167       IF (IALTH .EQ. 1) IALTH = 2
 168       IF (METH .EQ. MEO) GO TO 110
 169       CALL DCFODE (METH, ELCO, TESCO)
 170       MEO = METH
 171       IF (NQ .GT. MAXORD) GO TO 120
 172       IALTH = L
 173       IRET = 1
 174       GO TO 150
 175  110  IF (NQ .LE. MAXORD) GO TO 160
 176  120  NQ = MAXORD
 177       L = LMAX
 178       DO 125 I = 1,L
 179  125    EL(I) = ELCO(I,NQ)
 180       NQNYH = NQ*NYH
 181       RC = RC*EL(1)/EL0
 182       EL0 = EL(1)
 183       CONIT = 0.5D0/(NQ+2)
 184       DDN = DVNORM (N, SAVF, EWT)/TESCO(1,L)
 185       EXDN = 1.0D0/L
 186       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
 187       RH = MIN(RHDN,1.0D0)
 188       IREDO = 3
 189       IF (H .EQ. HOLD) GO TO 170
 190       RH = MIN(RH,ABS(H/HOLD))
 191       H = HOLD
 192       GO TO 175
 193 C-----------------------------------------------------------------------
 194 C DCFODE is called to get all the integration coefficients for the
 195 C current METH.  Then the EL vector and related constants are reset
 196 C whenever the order NQ is changed, or at the start of the problem.
 197 C-----------------------------------------------------------------------
 198  140  CALL DCFODE (METH, ELCO, TESCO)
 199  150  DO 155 I = 1,L
 200  155    EL(I) = ELCO(I,NQ)
 201       NQNYH = NQ*NYH
 202       RC = RC*EL(1)/EL0
 203       EL0 = EL(1)
 204       CONIT = 0.5D0/(NQ+2)
 205       GO TO (160, 170, 200), IRET
 206 C-----------------------------------------------------------------------
 207 C If H is being changed, the H ratio RH is checked against
 208 C RMAX, HMIN, and HMXI, and the YH array rescaled.  IALTH is set to
 209 C L = NQ + 1 to prevent a change of H for that many steps, unless
 210 C forced by a convergence or error test failure.
 211 C-----------------------------------------------------------------------
 212  160  IF (H .EQ. HOLD) GO TO 200
 213       RH = H/HOLD
 214       H = HOLD
 215       IREDO = 3
 216       GO TO 175
 217  170  RH = MAX(RH,HMIN/ABS(H))
 218  175  RH = MIN(RH,RMAX)
 219       RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
 220       R = 1.0D0
 221       DO 180 J = 2,L
 222         R = R*RH
 223         DO 180 I = 1,N
 224  180      YH(I,J) = YH(I,J)*R
 225       H = H*RH
 226       RC = RC*RH
 227       IALTH = L
 228       IF (IREDO .EQ. 0) GO TO 690
 229 C-----------------------------------------------------------------------
 230 C This section computes the predicted values by effectively
 231 C multiplying the YH array by the Pascal Triangle matrix.
 232 C RC is the ratio of new to old values of the coefficient  H*EL(1).
 233 C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
 234 C to force PJAC to be called, if a Jacobian is involved.
 235 C In any case, PJAC is called at least every MSBP steps.
 236 C-----------------------------------------------------------------------
 237  200  IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
 238       IF (NST .GE. NSLP+MSBP) IPUP = MITER
 239       TN = TN + H
 240       I1 = NQNYH + 1
 241       DO 215 JB = 1,NQ
 242         I1 = I1 - NYH
 243 Cdir$ ivdep
 244         DO 210 I = I1,NQNYH
 245  210      YH1(I) = YH1(I) + YH1(I+NYH)
 246  215    CONTINUE
 247 C-----------------------------------------------------------------------
 248 C Up to MAXCOR corrector iterations are taken.  A convergence test is
 249 C made on the R.M.S. norm of each correction, weighted by the error
 250 C weight vector EWT.  The sum of the corrections is accumulated in the
 251 C vector ACOR(i).  The YH array is not altered in the corrector loop.
 252 C-----------------------------------------------------------------------
 253  220  M = 0
 254       DO 230 I = 1,N
 255  230    Y(I) = YH(I,1)
 256       CALL F (NEQ, TN, Y, SAVF)
 257       NFE = NFE + 1
 258       IF (IPUP .LE. 0) GO TO 250
 259 C-----------------------------------------------------------------------
 260 C If indicated, the matrix P = I - h*el(1)*J is reevaluated and
 261 C preprocessed before starting the corrector iteration.  IPUP is set
 262 C to 0 as an indicator that this has been done.
 263 C-----------------------------------------------------------------------
 264       CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
 265       IPUP = 0
 266       RC = 1.0D0
 267       NSLP = NST
 268       CRATE = 0.7D0
 269       IF (IERPJ .NE. 0) GO TO 430
 270  250  DO 260 I = 1,N
 271  260    ACOR(I) = 0.0D0
 272  270  IF (MITER .NE. 0) GO TO 350
 273 C-----------------------------------------------------------------------
 274 C In the case of functional iteration, update Y directly from
 275 C the result of the last function evaluation.
 276 C-----------------------------------------------------------------------
 277       DO 290 I = 1,N
 278         SAVF(I) = H*SAVF(I) - YH(I,2)
 279  290    Y(I) = SAVF(I) - ACOR(I)
 280       DEL = DVNORM (N, Y, EWT)
 281       DO 300 I = 1,N
 282         Y(I) = YH(I,1) + EL(1)*SAVF(I)
 283  300    ACOR(I) = SAVF(I)
 284       GO TO 400
 285 C-----------------------------------------------------------------------
 286 C In the case of the chord method, compute the corrector error,
 287 C and solve the linear system with that as right-hand side and
 288 C P as coefficient matrix.
 289 C-----------------------------------------------------------------------
 290  350  DO 360 I = 1,N
 291  360    Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
 292       CALL SLVS (WM, IWM, Y, SAVF)
 293       IF (IERSL .LT. 0) GO TO 430
 294       IF (IERSL .GT. 0) GO TO 410
 295       DEL = DVNORM (N, Y, EWT)
 296       DO 380 I = 1,N
 297         ACOR(I) = ACOR(I) + Y(I)
 298  380    Y(I) = YH(I,1) + EL(1)*ACOR(I)
 299 C-----------------------------------------------------------------------
 300 C Test for convergence.  If M.gt.0, an estimate of the convergence
 301 C rate constant is stored in CRATE, and this is used in the test.
 302 C-----------------------------------------------------------------------
 303  400  IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
 304       DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
 305       IF (DCON .LE. 1.0D0) GO TO 450
 306       M = M + 1
 307       IF (M .EQ. MAXCOR) GO TO 410
 308       IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
 309       DELP = DEL
 310       CALL F (NEQ, TN, Y, SAVF)
 311       NFE = NFE + 1
 312       GO TO 270
 313 C-----------------------------------------------------------------------
 314 C The corrector iteration failed to converge.
 315 C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
 316 C the next try.  Otherwise the YH array is retracted to its values
 317 C before prediction, and H is reduced, if possible.  If H cannot be
 318 C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
 319 C-----------------------------------------------------------------------
 320  410  IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
 321       ICF = 1
 322       IPUP = MITER
 323       GO TO 220
 324  430  ICF = 2
 325       NCF = NCF + 1
 326       RMAX = 2.0D0
 327       TN = TOLD
 328       I1 = NQNYH + 1
 329       DO 445 JB = 1,NQ
 330         I1 = I1 - NYH
 331 Cdir$ ivdep
 332         DO 440 I = I1,NQNYH
 333  440      YH1(I) = YH1(I) - YH1(I+NYH)
 334  445    CONTINUE
 335       IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
 336       IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
 337       IF (NCF .EQ. MXNCF) GO TO 670
 338       RH = 0.25D0
 339       IPUP = MITER
 340       IREDO = 1
 341       GO TO 170
 342 C-----------------------------------------------------------------------
 343 C The corrector has converged.  JCUR is set to 0
 344 C to signal that the Jacobian involved may need updating later.
 345 C The local error test is made and control passes to statement 500
 346 C if it fails.
 347 C-----------------------------------------------------------------------
 348  450  JCUR = 0
 349       IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
 350       IF (M .GT. 0) DSM = DVNORM (N, ACOR, EWT)/TESCO(2,NQ)
 351       IF (DSM .GT. 1.0D0) GO TO 500
 352 C-----------------------------------------------------------------------
 353 C After a successful step, update the YH array.
 354 C Consider changing H if IALTH = 1.  Otherwise decrease IALTH by 1.
 355 C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
 356 C use in a possible order increase on the next step.
 357 C If a change in H is considered, an increase or decrease in order
 358 C by one is considered also.  A change in H is made only if it is by a
 359 C factor of at least 1.1.  If not, IALTH is set to 3 to prevent
 360 C testing for that many steps.
 361 C-----------------------------------------------------------------------
 362       KFLAG = 0
 363       IREDO = 0
 364       NST = NST + 1
 365       HU = H
 366       NQU = NQ
 367       DO 470 J = 1,L
 368         DO 470 I = 1,N
 369  470      YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
 370       IALTH = IALTH - 1
 371       IF (IALTH .EQ. 0) GO TO 520
 372       IF (IALTH .GT. 1) GO TO 700
 373       IF (L .EQ. LMAX) GO TO 700
 374       DO 490 I = 1,N
 375  490    YH(I,LMAX) = ACOR(I)
 376       GO TO 700
 377 C-----------------------------------------------------------------------
 378 C The error test failed.  KFLAG keeps track of multiple failures.
 379 C Restore TN and the YH array to their previous values, and prepare
 380 C to try the step again.  Compute the optimum step size for this or
 381 C one lower order.  After 2 or more failures, H is forced to decrease
 382 C by a factor of 0.2 or less.
 383 C-----------------------------------------------------------------------
 384  500  KFLAG = KFLAG - 1
 385       TN = TOLD
 386       I1 = NQNYH + 1
 387       DO 515 JB = 1,NQ
 388         I1 = I1 - NYH
 389 Cdir$ ivdep
 390         DO 510 I = I1,NQNYH
 391  510      YH1(I) = YH1(I) - YH1(I+NYH)
 392  515    CONTINUE
 393       RMAX = 2.0D0
 394       IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660
 395       IF (KFLAG .LE. -3) GO TO 640
 396       IREDO = 2
 397       RHUP = 0.0D0
 398       GO TO 540
 399 C-----------------------------------------------------------------------
 400 C Regardless of the success or failure of the step, factors
 401 C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
 402 C at order NQ - 1, order NQ, or order NQ + 1, respectively.
 403 C In the case of failure, RHUP = 0.0 to avoid an order increase.
 404 C The largest of these is determined and the new order chosen
 405 C accordingly.  If the order is to be increased, we compute one
 406 C additional scaled derivative.
 407 C-----------------------------------------------------------------------
 408  520  RHUP = 0.0D0
 409       IF (L .EQ. LMAX) GO TO 540
 410       DO 530 I = 1,N
 411  530    SAVF(I) = ACOR(I) - YH(I,LMAX)
 412       DUP = DVNORM (N, SAVF, EWT)/TESCO(3,NQ)
 413       EXUP = 1.0D0/(L+1)
 414       RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0)
 415  540  EXSM = 1.0D0/L
 416       RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
 417       RHDN = 0.0D0
 418       IF (NQ .EQ. 1) GO TO 560
 419       DDN = DVNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
 420       EXDN = 1.0D0/NQ
 421       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
 422  560  IF (RHSM .GE. RHUP) GO TO 570
 423       IF (RHUP .GT. RHDN) GO TO 590
 424       GO TO 580
 425  570  IF (RHSM .LT. RHDN) GO TO 580
 426       NEWQ = NQ
 427       RH = RHSM
 428       GO TO 620
 429  580  NEWQ = NQ - 1
 430       RH = RHDN
 431       IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
 432       GO TO 620
 433  590  NEWQ = L
 434       RH = RHUP
 435       IF (RH .LT. 1.1D0) GO TO 610
 436       R = EL(L)/L
 437       DO 600 I = 1,N
 438  600    YH(I,NEWQ+1) = ACOR(I)*R
 439       GO TO 630
 440  610  IALTH = 3
 441       GO TO 700
 442  620  IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
 443       IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
 444 C-----------------------------------------------------------------------
 445 C If there is a change of order, reset NQ, l, and the coefficients.
 446 C In any case H is reset according to RH and the YH array is rescaled.
 447 C Then exit from 690 if the step was OK, or redo the step otherwise.
 448 C-----------------------------------------------------------------------
 449       IF (NEWQ .EQ. NQ) GO TO 170
 450  630  NQ = NEWQ
 451       L = NQ + 1
 452       IRET = 2
 453       GO TO 150
 454 C-----------------------------------------------------------------------
 455 C Control reaches this section if 3 or more failures have occured.
 456 C If 10 failures have occurred, exit with KFLAG = -1.
 457 C It is assumed that the derivatives that have accumulated in the
 458 C YH array have errors of the wrong order.  Hence the first
 459 C derivative is recomputed, and the order is set to 1.  Then
 460 C H is reduced by a factor of 10, and the step is retried,
 461 C until it succeeds or H reaches HMIN.
 462 C-----------------------------------------------------------------------
 463  640  IF (KFLAG .EQ. -10) GO TO 660
 464       RH = 0.1D0
 465       RH = MAX(HMIN/ABS(H),RH)
 466       H = H*RH
 467       DO 645 I = 1,N
 468  645    Y(I) = YH(I,1)
 469       CALL F (NEQ, TN, Y, SAVF)
 470       NFE = NFE + 1
 471       DO 650 I = 1,N
 472  650    YH(I,2) = H*SAVF(I)
 473       IPUP = MITER
 474       IALTH = 5
 475       IF (NQ .EQ. 1) GO TO 200
 476       NQ = 1
 477       L = 2
 478       IRET = 3
 479       GO TO 150
 480 C-----------------------------------------------------------------------
 481 C All returns are made through this section.  H is saved in HOLD
 482 C to allow the caller to change H on the next step.
 483 C-----------------------------------------------------------------------
 484  660  KFLAG = -1
 485       GO TO 720
 486  670  KFLAG = -2
 487       GO TO 720
 488  680  KFLAG = -3
 489       GO TO 720
 490  690  RMAX = 10.0D0
 491  700  R = 1.0D0/TESCO(2,NQU)
 492       DO 710 I = 1,N
 493  710    ACOR(I) = ACOR(I)*R
 494  720  HOLD = H
 495       JSTART = 1
 496       RETURN
 497 C----------------------- END OF SUBROUTINE DSTODE ----------------------
 498       END