      SUBROUTINE COSET
C-----------------------------------------------------------------------
C COSET IS CALLED BY TSTEP AND SETS COEFFICIENTS FOR USE THERE.
C
C FOR EACH ORDER NQ, THE COEFFICIENTS IN EL ARE CALCULATED BY USE OF
C  THE GENERATING POLYNOMIAL LAMBDA(X), WITH COEFFICIENTS EL(I)
C      LAMBDA(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*(X**NQ).
C FOR THE BACKWARD DIFFERENTIATION FORMULAS,
C                    NQ
C      LAMBDA(X) = PRODUCT (1 + X/XI(I) ) .
C                   I = 1
C FOR THE ADAMS FORMULAS,
C                              NQ-1
C      (D/DX) LAMBDA(X) = C * PRODUCT (1 + X/XI(I) ) ,
C                              I = 1
C      LAMBDA(-1) = 0,    LAMBDA(0) = 1,
C WHERE C IS A NORMALIZATION CONSTANT.
C IN BOTH CASES, XI(I) IS DEFINED BY
C      H*XI(I) = T SUB N  -  T SUB (N-I)
C              = H + TAU(1) + TAU(2) + ... TAU(I-1).
C
C COSET ALSO SETS MAXDER, THE MAXIMUM ORDER OF THE FORMULAS
C AVAILABLE. CURRENTLY THIS IS 5 FOR THE BACKWARD DIFFERENTIATION
C FORMULAS, AND 12 FOR THE ADAMS FORMULAS.  TO USE DIFFERENT
C VALUES (.LE. 13),  CHANGE THE NUMBERS IN STATEMENTS 1 AND 2 BELOW.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH COSET USES THE FOLLOWING..
C   TAU    = A VECTOR OF LENGTH 13 CONTAINING THE PAST NQ VALUES
C            OF H.
C   EL     = A VECTOR OF LENGTH 13 IN WHICH COSET STORES THE
C            COEFFICIENTS FOR THE CORRECTOR FORMULA.
C   TQ     = A VECTOR OF LENGTH 5 IN WHICH COSET STORES CONSTANTS
C            USED FOR THE CONVERGENCE TEST, THE ERROR TEST, AND
C            SELECTION OF H AT A NEW ORDER.
C   LMAX   = MAXDER + 1, WHERE MAXDER IS THE MAXIMUM ORDER
C            AVAILABLE.  LMAX IS THE MAXIMUM NUMBER OF COLUMNS
C            OF THE Y ARRAY TO BE USED.
C   METH   = THE BASIC METHOD INDICATOR.
C   NQ     = THE CURRENT ORDER.
C   L      = NQ + 1, THE LENGTH OF THE VECTOR STORED IN EL, AND
C            THE NUMBER OF COLUMNS OF THE Y ARRAY BEING USED.
C   NQINDX = A COUNTER CONTROLLING THE FREQUENCY OF ORDER CHANGES.
C            AN ORDER CHANGE IS ABOUT TO BE CONSIDERED IF
C            NQINDX = 1.
C-----------------------------------------------------------------------
C CAUTION   NOT ALL MEMBERS OF EPCOM1 ARE USED IN THIS SUBROUTINE.
C-----------------------------------------------------------------------
      INTEGER I, IBACK, J, JP1, JSTART, KFLAG, L, LMAX, MAXDER,
     1        METH, MF,N, NQ, NQINDX, NQM1
      REAL EL, EPS, H, HMIN, SS, T, TAU, TQ, UROUND
      REAL AHDSS, CNQM1, CSUM, ELP, EM, EM0, FLOTI,
     1     FLOTL, FLOTNQ, HSUM, HSUM1, PROD, RXI, S, XI
      REAL CORTES
      REAL ONE, SIX, TWO, ZERO
      DIMENSION EM(13)
C
      COMMON /EPCOM1/ T,H,HMIN,HMAX,EPS,SS,UROUND,N,MF,KFLAG,JSTART
      COMMON /EPCM10/ TAU(13),EL(13),TQ(5),LMAX,METH,NQ,L,NQINDX
      DATA CORTES /0.1E0/
      DATA ONE  /1.0E0/, SIX /6.0E0/, TWO /2.0E0/, ZERO /0.0E0/
C
      AHDSS = ONE
      IF (SS .NE. ZERO) AHDSS = ABS(H)/SS
      FLOTL = FLOAT(L)
      NQM1 = NQ - 1
      GO TO (1, 2), METH
 1    MAXDER = 12
      GO TO 100
C
 2    MAXDER = 5
      GO TO 200
C
 100  IF (NQ .NE. 1) GO TO 110
      EL(1) = ONE
      EL(2) = ONE
      TQ(1) = ONE
      TQ(2) = TWO*AHDSS
      TQ(3) = SIX*TQ(2)
      TQ(5) = ONE
      GO TO 300
 110  HSUM = H
      EM(1) = ONE
      FLOTNQ = FLOTL - ONE
      DO 115 I = 2,L
 115    EM(I) = ZERO
      DO 150 J = 1,NQM1
        IF ((J .NE. NQM1) .OR. (NQINDX .NE. 1)) GO TO 130
        S = ONE
        CSUM = ZERO
        DO 120 I = 1,NQM1
          CSUM = CSUM + S*EM(I)/FLOAT(I+1)
 120      S = -S
        TQ(1) = AHDSS*EM(NQM1)/(FLOTNQ*CSUM)
 130    RXI = H/HSUM
        DO 140 IBACK = 1,J
          I = (J + 2) - IBACK
 140      EM(I) = EM(I) + EM(I-1)*RXI
 150    HSUM = HSUM + TAU(J)
C COMPUTE INTEGRAL FROM -1 TO 0 OF POLYNOMIAL AND OF X TIMES IT. -------
      S = ONE
      EM0 = ZERO
      CSUM = ZERO
      DO 160 I = 1,NQ
        FLOTI = FLOAT(I)
        EM0 = EM0 + S*EM(I)/FLOTI
        CSUM = CSUM + S*EM(I)/(FLOTI+1)
 160    S = -S
C IN EL, FORM COEFFICIENTS OF NORMALIZED INTEGRATED POLYNOMIAL. --------
      S = ONE/EM0
      EL(1) = ONE
      DO 170 I = 1,NQ
  170   EL(I+1) = S*EM(I)/FLOAT(I)
      XI = HSUM/H
      TQ(2) = AHDSS*XI*EM0/CSUM
      TQ(5) = XI/EL(L)
      IF (NQINDX .NE. 1) GO TO 300
C FOR HIGHER ORDER CONTROL CONSTANT, MULTIPLY POLYNOMIAL BY 1+X/XI(Q). -
      RXI = ONE/XI
      DO 180 IBACK = 1,NQ
        I = (L + 1) - IBACK
 180    EM(I) = EM(I) + EM(I-1)*RXI
C COMPUTE INTEGRAL OF POLYNOMIAL. --------------------------------------
      S = ONE
      CSUM = ZERO
      DO 190 I = 1,L
      CSUM = CSUM + S*EM(I)/FLOAT(I+1)
 190    S = -S
      TQ(3) = AHDSS*FLOTL*EM0/CSUM
      GO TO 300
C
 200  DO 210 I = 3,L
 210    EL(I) = ZERO
      EL(1) = ONE
      EL(2) = ONE
      HSUM = H
      HSUM1 = ZERO
      PROD = ONE
      RXI = ONE
      IF (NQ .EQ. 1) GO TO 240
      DO 230 J = 1,NQM1
C IN EL, CONSTRUCT COEFFICIENTS OF (1+X/XI(1))*...*(1+X/XI(J+1)). ------
        HSUM = HSUM + TAU(J)
        HSUM1 = HSUM1 + TAU(J)
        PROD = PROD*(HSUM/HSUM1)
        RXI = H/HSUM
        JP1 = J + 1
        DO 220 IBACK = 1,JP1
          I = (J + 3) - IBACK
 220      EL(I) = EL(I) + EL(I-1)*RXI
 230    CONTINUE
 240  TQ(2) = AHDSS*EL(2)*(ONE + PROD)
      TQ(5) = (ONE + PROD)/EL(L)
      IF (NQINDX .NE. 1) GO TO 300
      CNQM1 = RXI/EL(L)
      ELP = EL(2) - RXI
      TQ(1) = AHDSS*ELP/CNQM1
      HSUM = HSUM + TAU(NQ)
      RXI = H/HSUM
      ELP = EL(2) + RXI
      TQ(3) = AHDSS*ELP*RXI*(ONE + PROD)*(FLOTL + ONE)
 300  TQ(4) = CORTES*TQ(2)
      LMAX = MAXDER + 1
      RETURN
C-----------------------  END OF SUBROUTINE COSET  ---------------------
      END
