      SUBROUTINE TSTEP(Y,N0,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,
     1                   DIFFUN,PEDERV)
C-----------------------------------------------------------------------
C TSTEP PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
C COMMUNICATION WITH TSTEP IS VIA THE FOLLOWING VARIABLES..
C
C   Y       AN N0 BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
C             AND THEIR SCALED DERIVATIVES.  LMAX IS CURRENTLY 6 FOR
C             THE VARIABLE STEP BACKWARD DIFFERENTIATION FORMULAS,
C             AND 13 FOR THE VARIABLE STEP ADAMS FORMULAS.
C             (LMAX -1) = MAXDER, THE MAXIMUM ORDER USED.
C             SEE SUBROUTINE COSET.  Y(I,J+1) CONTAINS THE
C             J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J)
C             FOR J = 0,1,...,NQ, WHERE NQ IS THE CURRENT ORDER.
C   N0      A CONSTANT INTEGER .GE. N, USED FOR DIMENSIONING
C             PURPOSES.
C   T       THE INDEPENDENT VARIABLE, UPDATED ON EACH STEP TAKEN.
C   H       THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
C             H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING
C             THE SOLUTION OF THE PROBLEM. H CAN BE EITHER POSITIVE
C             OR NEGATIVE, BUT ITS SIGN MUST REMAIN CONSTANT
C             THROUGHOUT THE PROBLEM RUN.
C   HMIN,   THE MINIMUM AND MAXIMUM ABSOLUTE VALUES OF THE STEP
C    HMAX     SIZE TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT
C             ANY TIME, BUT THE CHANGE WILL NOT TAKE EFFECT UNTIL THE
C             NEXT CHANGE IN H IS MADE.
C   EPS     THE RELATIVE ERROR BOUND.  SEE DESCRIPTION IN
C             SUBROUTINE DRIVE.
C   SS      THE SIZE OF THE TIME INTERVAL TO BE USED FOR ERROR
C             CONTROL.  A DEFAULT VALUE OF 0 IS USED TO PRODUCE
C             CONTROL OF ERROR PER STEP.  SEE SUBROUTINE DRIVE.
C   UROUND  THE UNIT OF ROUNDOFF FOR THE COMPUTER BEING USED.
C   N       THE NUMBER OF FIRST ORDER ORDINARY DIFFERENTIAL
C             EQUATIONS BEING SOLVED.
C   MF      THE METHOD FLAG.  SEE DESCRIPTION IN SUBROUTINE DRIVE.
C   KFLAG   A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
C                     0  THE STEP WAS SUCCESFUL.
C                    -1  THE REQUESTED ERROR COULD NOT BE ACHIEVED
C                          WITH ABS(H) = HMIN.
C                    -2  THE REQUESTED ERROR IS SMALLER THAN CAN
C                          BE HANDLED FOR THIS PROBLEM.
C                    -3  CORRECTOR CONVERGENCE COULD NOT BE
C                          ACHIEVED FOR ABS(H) = HMIN.
C             ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND
C             THE Y ARRAY ARE AS OF THE BEGINNING OF THE LAST
C             STEP AND H IS THE LAST STEP SIZE ATTEMPTED.
C   JSTART  AN INTEGER USED ON INPUT AND OUTPUT.
C             ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS..
C                     0  PERFORM THE FIRST STEP.
C                 .GT.0  TAKE A NEW STEP CONTINUING FROM THE LAST.
C                 .LT.0  TAKE THE NEXT STEP WITH A NEW VALUE OF
C                          H AND/OR MF.
C             ON EXIT, JSTART IS SET TO NQ, THE CURRENT ORDER OF THE
C             METHOD.
C   YMAX      AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL
C               ERRORS IN Y ARE COMPARED.
C   ERROR     AN ARRAY OF N ELEMENTS.  ERROR(I)/TQ(2) IS THE
C               ESTIMATED LOCAL ERROR IN Y(I) PER SS UNITS OF
C               T OR PER STEP (OF SIZE H).
C   SAVE1,    TWO ARRAYS FOR WORKING STORAGE,
C     SAVE2     EACH OF LENGTH N.
C   PW        A BLOCK OF LOCATIONS USED FOR THE PARTIAL DERIVATIVES
C               OF F WITH RESPECT TO Y, IF MITER IS NOT O.  SEE
C               DESCRIPTION IN SUBROUTINE DRIVE.
C   IPIV      AN INTEGER ARRAY OF LENGTH N, WHICH IS USED FOR PIVOT
C               INFORMATION FOR THE LINEAR ALGEBRAIC SYSTEM IN THE
C               CORRECTION PROCESS, WHEN MITER = 1 OR 2.
C
C THE COMMON BLOCK EPCM10, DECLARED BELOW, IS PRIMARILY INTENDED
C FOR INTERNAL USE, BUT IT CAN BE ACCESSED EXTERNALLY.
C-----------------------------------------------------------------------
      INTEGER N0
      INTEGER IPIV, JSTART, KFLAG, L, LMAX, METH, MF, N, NFE, NJE,
     1        NQ, NQINDX, NQUSED, NSQ, NSTEP
      INTEGER I, IBACK, IER, IREDO, J, J1, J2, M, MFOLD, MIO,
     1        MITER, MITER1, NEWJ, NSTEPJ
      INTEGER ISTEPJ, KFC, KFH, MAXCOR
      REAL Y
      REAL EL, EPS, EPSJ, ERROR, H, HMAX, HMIN, HUSED, PW,
     1     SAVE1, SAVE2, SS, T, TAU, TQ, UROUND, YMAX
      REAL BND, CNQUOT, CON, CONP, CRATE, D, DRC,
     1     D1, E, EDN, ETA, ETAMAX, ETAMIN, ETAQ, ETAQM1,
     2     ETAQP1, EUP, FLOTL, FLOTN, HOLD, HRL1, PHRL1,
     3     PRL1, R, RC, RL1, R0, R1, TOLD
      REAL ADDON, BIAS1, BIAS2, BIAS3, CRDOWN, DELRC,
     1     ETACF, ETAMXF, ETAMX1, ETAMX2,
     2     ETAMX3, ONEPSM, SHORT, THRESH
      REAL ONE, PT5, ZERO
      DIMENSION Y(N0,13)
      DIMENSION YMAX(1),ERROR(1),SAVE1(1),SAVE2(1),PW(1),IPIV(1)
      EXTERNAL DIFFUN,PEDERV
C
      COMMON /XXXSTP/ CONP,CRATE,D1,ETA,ETAMAX,HOLD,HRL1,IREDO,MFOLD,
     1                MITER,MITER1,NEWJ,NSTEPJ,PRL1,RC
      COMMON /EPCOM1/ T,H,HMIN,HMAX,EPS,SS,UROUND,N,MF,KFLAG,JSTART
      COMMON /EPCOM8/ EPSJ,NSQ
      COMMON /EPCOM9/ HUSED,NQUSED,NSTEP,NFE,NJE
      COMMON /EPCM10/ TAU(13),EL(13),TQ(5),LMAX,METH,NQ,L,NQINDX
C
      DATA ISTEPJ /20/, KFC /-3/, KFH /-7/, MAXCOR /3/
      DATA ADDON  /1.0E-6/,    BIAS1  /2.5E1/, BIAS2  /2.5E1/,
     1     BIAS3  /1.0E2/,     CRDOWN /0.1E0/, DELRC  /0.3E0/,
     2     ETACF  /0.25E0/,    ETAMIN /0.1E0/, ETAMXF /0.2E0/,
     3     ETAMX1 /1.0E4/,     ETAMX2 /1.0E1/, ETAMX3 /1.5E0/,
     4     ONEPSM /1.00001E0/, SHORT  /0.1E0/, THRESH /1.3E0/
      DATA ONE /1.0E0/, PT5 /0.5E0/, ZERO /0.0E0/
C
      KFLAG = 0
      TOLD = T
      FLOTN = FLOAT(N)
      IF (JSTART .GT. 0) GO TO 200
      IF (JSTART .NE. 0) GO TO 150
C-----------------------------------------------------------------------
C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL
C DERIVATIVES ARE CALCULATED.  ETAMAX IS THE MAXIMUM RATIO BY
C WHICH H CAN BE INCREASED IN A SINGLE STEP.  IT IS 1.E04 FOR THE
C FIRST STEP TO COMPENSATE FOR THE SMALL INITIAL H, THEN 10 FOR
C THE NEXT 10 STEPS, AND THEN 1.5 THEREAFTER.  IF A FAILURE
C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), ETAMAX IS SET AT 1
C FOR THE NEXT INCREASE.  ETAMIN = .1 IS THE MINIMUM RATIO BY WHICH
C H CAN BE REDUCED ON ANY RETRY OF A STEP.
C-----------------------------------------------------------------------
      CALL DIFFUN (N, T, Y, SAVE1)
      DO 110 I = 1,N
 110    Y(I,2) = H*SAVE1(I)
      METH = MF/10
      MITER = MF - 10*METH
      MITER1 = MITER + 1
      MFOLD = MF
      NQ = 1
      L = 2
      TAU(1) = H
      PRL1 = ONE
      RC = ZERO
      ETAMAX = ETAMX1
      NQINDX = 2
      NSTEP = 0
      NSTEPJ = 0
      NFE = 1
      NJE = 0
      GO TO 200
C-----------------------------------------------------------------------
C IF THE USER HAS CHANGED H, THEN Y MUST BE RESCALED.  IF THE
C USER HAS CHANGED MITER, THEN NEWJ IS SET TO MITER TO FORCE
C THE PARTIAL DERIVATIVEES TO BE UPDATED, IF THEY ARE BEING USED.
C-----------------------------------------------------------------------
 150  IF (MF .EQ. MFOLD) GO TO 170
      MIO = MITER
      METH = MF/10
      MITER = MF - 10*METH
      MFOLD = MF
      IF (MITER .EQ. MIO) GO TO 170
      NEWJ = MITER
      MITER1 = MITER + 1
 170  IF (H .EQ. HOLD) GO TO 200
      ETA = H/HOLD
      H = HOLD
      IREDO = 3
      GO TO 185
  180 ETA = AMAX1(ETA,HMIN/ABS(H),ETAMIN)
  185 ETA = AMIN1(ETA,HMAX/ABS(H),ETAMAX)
      R1 = ONE
      DO 190 J = 2,L
        R1 = R1*ETA
        DO 190 I = 1,N
 190      Y(I,J) = Y(I,J)*R1
      H = H*ETA
      RC = RC*ETA
      IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
C MULTIPLYING THE Y ARRAY BY THE PASCAL TRIANGLE MATRIX.  THEN
C COSET IS CALLED TO OBTAIN EL, THE VECTOR OF COEFFICIENTS OF
C LENGTH NQ + 1.  RC IS THE RATIO OF NEW TO OLD VALUES OF THE
C COEFFICIENT H/EL(2).  WHEN RC DIFFERS FROM 1 BY MORE THAN
C DELRC, NEWJ IS SET TO MITER TO FORCE THE PARTIAL DERIVATIVES
C TO BE UPDATED, IF USED.  DELRC IS 0.3.  IN ANY CASE, THE PARTIAL
C DERIVATIVES ARE UPDATED AT LEAST EVERY 20-TH STEP.
C-----------------------------------------------------------------------
 200  T = T + H
      DO 210 J1 = 1,NQ
        DO 210 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 210 I = 1,N
 210        Y(I,J) = Y(I,J) + Y(I,J+1)
      CALL COSET
      BND = FLOTN*(TQ(4)*EPS)**2
      RL1 = ONE/EL(2)
      RC = RC*(RL1/PRL1)
      PRL1 = RL1
      IF (NSTEP .GE. NSTEPJ+ISTEPJ) NEWJ = MITER
      DRC = ABS(RC-ONE)
      IF (DRC .LE. DELRC) GO TO 215
      NEWJ = MITER
      CRATE = ONE
      RC = ONE
      GO TO 220
  215 IF ((MITER .NE. 0) .AND. (DRC .NE. ZERO)) CRATE = ONE
C-----------------------------------------------------------------------
C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN.  A CONVERGENCE TEST IS MADE
C ON THE ROOT MEAN SQUARE NORM OF EACH CORRECTION, USING BND, WHICH
C IS DEPENDENT ON EPS.  THE SUM OF THE CORRECTIONS IS ACCUMULATED IN
C THE VECTOR ERROR.  THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR
C LOOP.  THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1.
C-----------------------------------------------------------------------
 220  DO 230 I = 1,N
 230    ERROR(I) = ZERO
      M = 0
      CALL DIFFUN (N, T, Y, SAVE2)
      NFE = NFE + 1
      IF (NEWJ .LE. 0) GO TO 290
C-----------------------------------------------------------------------
C IF INDICATED, THE MATRIX P = I - H*RL1*J IS REEVALUATED BEFORE
C STARTING THE CORRECTOR ITERATION.  NEWJ IS SET TO 0 AS AN
C INDICATOR THAT THIS HAS BEEN DONE.  IF MITER = 1 OR 2, P IS
C COMPUTED AND PROCESSED IN PSET.  IF MITER = 3, THE MATRIX  IS
C P = I - H*RL1*D, WHERE D IS A DIAGONAL MATRIX.  RL1 IS 1/EL(2).
C-----------------------------------------------------------------------
      NEWJ = 0
      RC = ONE
      NJE = NJE + 1
      NSTEPJ = NSTEP
      GO TO (250, 240, 260), MITER
 240  NFE = NFE + N
 250  CON = -H*RL1
      CALL PSET(Y,N0,CON,MITER,IER,YMAX,SAVE1,SAVE2,PW,IPIV,
     1            DIFFUN,PEDERV)
      IF (IER .NE. 0) GO TO 420
      GO TO 350
 260  R = RL1*SHORT
      DO 270 I = 1,N
 270    PW(I) = Y(I,1) + R*(H*SAVE2(I) - Y(I,2))
      CALL DIFFUN(N, T, PW, SAVE1)
      NFE = NFE + 1
      HRL1 = H*RL1
      DO 280 I = 1,N
        R0 = H*SAVE2(I) - Y(I,2)
        PW(I) = ONE
        D = SHORT*R0 - H*(SAVE1(I) - SAVE2(I))
        SAVE1(I) = ZERO
        IF (ABS(R0) .LT. UROUND*YMAX(I)) GO TO 280
        IF (ABS(D) .EQ. ZERO) GO TO 420
        PW(I) = SHORT*R0/D
        SAVE1(I) = PW(I)*RL1*R0
 280    CONTINUE
      GO TO 370
 290  GO TO (295, 350, 350, 310), MITER1
C-----------------------------------------------------------------------
C IN THE CASE OF FUNCTIONAL ITERATION, Y IS UPDATED DIRECTLY FROM
C THE RESULT OF THE LAST DIFFUN CALL.
C-----------------------------------------------------------------------
 295  D = ZERO
      DO 300 I = 1,N
        R = RL1*(H*SAVE2(I) - Y(I,2))
        D = D + ((R - ERROR(I))/YMAX(I))**2
        SAVE1(I) = Y(I,1) + R
 300    ERROR(I) = R
      GO TO 400
C-----------------------------------------------------------------------
C IN THE CASE OF A CHORD METHOD, THE RESIDUAL -G(Y SUB N(M))
C IS COMPUTED AND THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE
C AND P AS COEFFICIENT MATRIX IS SOLVED, USING THE LU DECOMPOSITION
C OF P IF MITER = 1 OR 2.  IF MITER = 3 THE SCALAR H*RL1 IS UPDATED.
C-----------------------------------------------------------------------
 310  PHRL1 = HRL1
      HRL1 = H*RL1
      IF (HRL1 .EQ. PHRL1) GO TO 330
      R = HRL1/PHRL1
      DO 320 I = 1,N
        D = ONE - R*(ONE - ONE/PW(I))
        IF (ABS(D) .EQ. ZERO) GO TO 440
 320    PW(I) = ONE/D
 330  DO 340 I = 1,N
 340    SAVE1(I) = PW(I)*(RL1*H*SAVE2(I) - (RL1*Y(I,2) + ERROR(I)))
      GO TO 370
 350  DO 360 I = 1,N
 360    SAVE1(I) = RL1*H*SAVE2(I) - (RL1*Y(I,2) + ERROR(I))
      CALL SOL (N, N0, PW, SAVE1, IPIV)
 370  D = ZERO
      DO 380 I = 1,N
        ERROR(I) = ERROR(I) + SAVE1(I)
        D = D + (SAVE1(I)/YMAX(I))**2
 380    SAVE1(I) = Y(I,1) + ERROR(I)
C-----------------------------------------------------------------------
C TEST FOR CONVERGENCE.  IF M .GT. 0, AN ESTIMATE OF THE SQUARE OF
C THE CONVERGENCE RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED
C IN THE TEST.
C-----------------------------------------------------------------------
  400 IF (M .NE. 0) CRATE = AMAX1(CRDOWN*CRATE,D/D1)
      IF (D*AMIN1(ONE,CRATE) .LE. BND) GO TO 450
      D1 = D
      M = M + 1
      IF (M .EQ. MAXCOR) GO TO 410
      CALL DIFFUN (N, T, SAVE1, SAVE2)
      GO TO (295, 350, 350, 310), MITER1
C-----------------------------------------------------------------------
C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. IF PARTIAL
C DERIVATIVES ARE INVOLVED BUT ARE NOT UP TO DATE, THEY ARE
C REEVALUATED FOR THE NEXT TRY.  OTHERWISE THE Y ARRAY IS RESTORED
C TO ITS VALUES BEFORE PREDICTION, AND H IS REDUCED,
C IF POSSIBLE.  IF NOT, A NO-CONVERGENCE EXIT IS TAKEN.
C-----------------------------------------------------------------------
 410  NFE = NFE + MAXCOR - 1
      IF (NEWJ .EQ. -1) GO TO 440
 420  T = TOLD
      ETAMAX = ONE
      DO 430 J1 = 1,NQ
        DO 430 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 430 I = 1,N
 430        Y(I,J) = Y(I,J) - Y(I,J+1)
      IF (ABS(H) .LE. HMIN*ONEPSM) GO TO 680
      ETA = ETACF
      IREDO = 1
      GO TO 180
 440  NEWJ = MITER
      GO TO 220
C-----------------------------------------------------------------------
C THE CORRECTOR HAS CONVERGED.  NEWJ IS SET TO -1 IF PARTIAL
C DERIVATIVES WERE USED, TO SIGNAL THAT THEY MAY NEED UPDATING ON
C SUBSEQUENT STEPS.  THE ERROR TEST IS MADE AND CONTROL PASSES TO
C STATEMENT 500 IF IT FAILS.
C-----------------------------------------------------------------------
 450  IF (MITER .NE. 0) NEWJ = -1
      NFE = NFE + M
      D = ZERO
      DO 460 I = 1,N
 460    D = D + (ERROR(I)/YMAX(I))**2
      E = FLOTN*(TQ(2)*EPS)**2
      IF (D .GT. E) GO TO 500
C-----------------------------------------------------------------------
C AFTER A SUCCESSFUL STEP, THE Y ARRAY, TAU, NSTEP, AND NQINDX ARE
C UPDATED, AND A NEW VALUE OF H AT ORDER NQ IS COMPUTED.
C THE VECTOR TAU CONTAINS THE NQ + 1 MOST RECENT VALUES OF H.
C A CHANGE IN NQ UP OR DOWN BY 1 IS CONSIDERED IF NQINDX = 0.
C IF NQINDX = 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED
C FOR USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
C A CHANGE IN H OR NQ IS MADE ONLY OF THE INCREASE IN H
C IS BY A FACTOR OF AT LEAST 1.3.
C IF NOT, NQINDX IS SET TO 2 TO PREVENT TESTING FOR THAT MANY
C STEPS.  IF NQ IS CHANGED, NQINDX IS SET TO NQ + 1 (NEW VALUE).
C-----------------------------------------------------------------------
      KFLAG = 0
      IREDO = 0
      NSTEP = NSTEP + 1
      HUSED = H
      NQUSED = NQ
      DO 470 IBACK = 1,NQ
        I = L - IBACK
 470    TAU(I+1) = TAU(I)
      TAU(1) = H
      DO 480 J = 1,L
        DO 480 I = 1,N
 480      Y(I,J) = Y(I,J) + ERROR(I)*EL(J)
      NQINDX = NQINDX - 1
      IF ((L .EQ. LMAX) .OR. (NQINDX .NE. 1)) GO TO 495
      DO 490 I = 1,N
 490    Y(I,LMAX) = ERROR(I)
      CONP = TQ(5)
 495  IF (ETAMAX .NE. ONE) GO TO 520
      IF (NQINDX .LT. 2) NQINDX = 2
      GO TO 690
C-----------------------------------------------------------------------
C THE ERROR TEST FAILED.  KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
C T AND THE Y ARRAY ARE RESTORED TO THEIR PREVIOUS VALUES.  A NEW
C H FOR A RETRY OF THE STEP IS COMPUTED.  THE ORDER IS KEPT FIXED.
C-----------------------------------------------------------------------
 500  KFLAG = KFLAG - 1
      T = TOLD
      DO 510 J1 = 1,NQ
        DO 510 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 510 I = 1,N
 510        Y(I,J) = Y(I,J) - Y(I,J+1)
      NEWJ = MITER
      ETAMAX = ONE
      IF (ABS(H) .LE. HMIN*ONEPSM) GO TO 660
      IF (KFLAG .LE. KFC) GO TO 630
      IREDO = 2
C COMPUTE RATIO OF NEW H TO CURRENT H AT THE CURRENT ORDER. ------------
  520 FLOTL = FLOAT(L)
      ETAQ = ONE/((BIAS2*D/E)**(PT5/FLOTL) + ADDON)
      IF ((NQINDX .NE. 0) .OR. (KFLAG .NE. 0)) GO TO 580
      ETAQM1 = ZERO
      IF (NQ .EQ. 1) GO TO 540
C COMPUTE RATIO OF NEW H TO CURRENT H AT THE CURRENT ORDER LESS ONE. ---
      D = ZERO
      DO 530 I = 1,N
 530    D = D + (Y(I,L)/YMAX(I))**2
      EDN = FLOTN*(TQ(1)*EPS)**2
      ETAQM1 = ONE/((BIAS1*D/EDN)**(PT5/(FLOTL - ONE)) + ADDON)
 540  ETAQP1 = ZERO
      IF (L .EQ. LMAX) GO TO 560
C COMPUTE RATIO OF NEW H TO CURRENT H AT CURRENT ORDER PLUS ONE. -------
      CNQUOT = (TQ(5)/CONP)*(H/TAU(2))**L
      D = ZERO
      DO 550 I = 1,N
 550    D = D + ((ERROR(I) - CNQUOT*Y(I,LMAX))/YMAX(I))**2
      EUP = FLOTN*(TQ(3)*EPS)**2
      ETAQP1 = ONE/((BIAS3*D/EUP)**(PT5/(FLOTL + ONE)) + ADDON)
 560  NQINDX = 2
      IF (ETAQ .GE. ETAQP1) GO TO 570
      IF (ETAQP1 .GT. ETAQM1) GO TO 600
      GO TO 590
 570  IF (ETAQ .LT. ETAQM1) GO TO 590
 580  IF ((ETAQ .LT. THRESH) .AND. (KFLAG .EQ. 0)) GO TO 690
      ETA = ETAQ
      IF ((KFLAG .LE. -2) .AND. (ETA .GT. ETAMXF)) ETA = ETAMXF
      GO TO 180
 590  IF (ETAQM1 .LT. THRESH) GO TO 690
      CALL ADJUST (Y, N0)
      L = NQ
      NQ = NQ - 1
      ETA = ETAQM1
      NQINDX = L
      GO TO 180
 600  IF (ETAQP1 .LT. THRESH) GO TO 690
      NQ = L
      ETA = ETAQP1
      L = L + 1
      DO 610 I = 1,N
 610    Y(I,L) = ZERO
      NQINDX = L
      GO TO 180
C-----------------------------------------------------------------------
C CONTROL REACHES THIS SECTION IF 3 OR MORE CONSECUTIVE FAILURES
C HAVE OCCURRED.  IT IS ASSUMED THAT THE ELEMENTS OF THE Y ARRAY
C HAVE ACCUMULATED ERRORS OF THE WRONG ORDER.  THE ORDER IS REDUCED
C BY ONE, IF POSSIBLE.  THEN H IS REDUCED BY A FACTOR OF 0.1 AND
C THE STEP IS RETRIED.  AFTER A TOTAL OF 7 CONSECUTIVE FAILURES,
C AN EXIT IS TAKEN WITH KFLAG = -2.
C-----------------------------------------------------------------------
 630  IF (KFLAG .EQ. KFH) GO TO 670
      IF (NQ .EQ. 1) GO TO 640
      ETA = ETAMIN
      CALL ADJUST (Y, N0)
      L = NQ
      NQ = NQ - 1
      NQINDX = L
      GO TO 180
  640 ETA = AMAX1(ETAMIN,HMIN/ABS(H))
      ETA = AMIN1(ETA,HMAX/ABS(H))
      H = H*ETA
      CALL DIFFUN (N, T, Y, SAVE1)
      NFE = NFE + 1
      DO 650 I = 1,N
 650    Y(I,2) = H*SAVE1(I)
      NQINDX = 10
      GO TO 200
C-----------------------------------------------------------------------
C ALL RETURNS ARE MADE THROUGH THIS SECTION.  H IS SAVED IN HOLD
C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
C-----------------------------------------------------------------------
 660  KFLAG = -1
      GO TO 700
 670  KFLAG = -2
      GO TO 700
 680  KFLAG = -3
      GO TO 700
 690  ETAMAX = ETAMX3
      IF (NSTEP .LE. 10) ETAMAX = ETAMX2
 700  HOLD = H
      JSTART = NQ
      RETURN
C-----------------------  END OF SUBROUTINE TSTEP  ---------------------
      END
