      SUBROUTINE CHEBY1 (A, B, F, G, PHI, EPS, ITNO, MXITER, L, M,
     1                   P, Q, ERROR, IERR, LP1, MP1, LPM, N, NP1,
     2                   X, XVAL, C, D, ERR, H)
C     ------------------------------------------------------------------
      DOUBLE PRECISION A, B, F, G, PHI, EPS, ERROR
      DOUBLE PRECISION P(LP1), Q(MP1), X(NP1), XVAL(NP1), C(NP1,NP1),
     1          D(NP1), ERR(NP1), H(NP1)
      DOUBLE PRECISION B1, C0, DEL, DN, DNP1, EPS0, HALF, H1,
     1          OLDERR, ONE, PI, SIGN, SUM, TAU, TEMPL, TEN, TEST,
     2          U, XI, XLB, XM1, Y, Y2, Y3, Z, ZERO, ZZ, Z1, Z2, Z3
      EXTERNAL F, G, PHI
C     -------------------
      DATA PI/3.14159265358979323846264338328D0/
      DATA ZERO/0.D0/, HALF/.5D0/, ONE/1.D0/, TEN/10.D0/
      DATA EPS0/1.D-2/, TAU/.015D0/, C0/.0625D0/
C     -------------------
      ERROR = ZERO
      IF (EPS .LE. ZERO .OR. EPS .GE. EPS0) GO TO 200
      IERR = 0
C
      ITNO = 1
      XLB = ZERO
      DN = N
      DNP1 = NP1
C
      DO 10 I = 1,LP1
   10 P(I) = ZERO
      DO 11 I = 1,MP1
   11 Q(I) = ZERO
      Q(1) = ONE
C
C     COMPUTE INITIAL APPROXIMATIONS OF THE CRITICAL POINTS
C
      X(1) = A
      X(NP1) = B
      K = N/2
      IF (K .LE. 0) GO TO 30
      B1 = HALF*(B - A)
      XM1 = HALF*(A + B)
      DO 20 I = 1,K
      XI = I
      Z = -B1*DCOS(PI*(XI/DN))
      X(I+1) = Z + XM1
      II = NP1 - I
   20 X(II) = XM1 - Z
C
C     EVALUATE PHI AT THE CRITICAL POINTS
C
   30 DO 31 I = 1,NP1
   31 XVAL(I) = PHI(X(I))
      KOUNT = 1
C
C     SET UP THE LINEAR EQUATIONS
C
   40 K = L + 2
      SIGN = ONE
      DO 45 I = 1,NP1
      SIGN = -SIGN
      C(I,1) = ONE
      IF (L .LE. 0) GO TO 42
      DO 41 J = 2,LP1
   41 C(I,J) = C(I,J-1)*XVAL(I)
   42 D(I) = F(X(I))
      IF (M .LE. 0) GO TO 44
      TEMPL = SIGN*XLB*G(X(I)) - D(I)
      C(I,K) = XVAL(I)*TEMPL
      IF (K .GT. LPM) GO TO 44
      DO 43 J = K,LPM
   43 C(I,J+1) = C(I,J)*XVAL(I)
   44 C(I,NP1) = SIGN*G(X(I))
   45 CONTINUE
C
C     SOLVE THE EQUATIONS CX = D AND STORE THE RESULTS IN D
C
      CALL DPSLV (NP1, 1, C, NP1, D, NP1, IERR)
      IF (IERR .NE. 0) GO TO 220
      IF (KOUNT .GT. 1) GO TO 50
C
C     REDEFINE THE EQUATIONS AND SOLVE
C
      XLB = (D(NP1) + XLB*DN)/DNP1
      IF (M .LE. 0) GO TO 61
      KOUNT = 2
      GO TO 40
C
   50 TEST = DABS(XLB - D(NP1))
      XLB = (D(NP1) + XLB*DN)/DNP1
      KOUNT = KOUNT + 1
      IF (KOUNT .LE. 4 .AND. TEST .GT. EPS0*DABS(XLB)) GO TO 40
C
C     STORE THE RESULTS IN P AND Q
C
      DO 60 I = 2,MP1
      LPI = L + I
   60 Q(I) = D(LPI)
   61 DO 62 I = 1,LP1
   62 P(I) = D(I)
C
C     SEARCH FOR NEW CRITICAL POINTS
C
      OLDERR = ERROR
      ERROR = ZERO
      Z1 = ZERO
      U = ONE
      IF (XLB .LT. ZERO) U = -U
C
      IF (N .GT. 1) GO TO 70
      H(1) = TAU*(X(2) - X(1))
      H(2) = -H(1)
      GO TO 72
   70 DO 71 I = 2,N
   71 H(I) = TAU*(X(I+1) - X(I-1))
      H(1) = HALF*H(2)
      H(NP1) = -HALF*H(N)
   72 CONTINUE
C
      DO 92 I = 1,NP1
      Y2 = X(I)
      H1 = H(I)
      Y3 = Y2 + H1
      CALL CERR(Y2, F(Y2), G(Y2), PHI(Y2), DEL, IERR, L, LP1, M, NP1, D)
      IF (IERR .NE. 0) RETURN
      Z2 = U*DEL
      CALL CERR(Y3, F(Y3), G(Y3), PHI(Y3), DEL, IERR, L, LP1, M, NP1, D)
      IF (IERR .NE. 0) RETURN
      Z3 = U*DEL
      IF (Z2 .LT. Z3) GO TO 80
      H1 = -H1
      Z = Z3
      Z3 = Z2
      Z2 = Z
      Y = Y3
      Y3 = Y2
      Y2 = Y
C
   80 Y = Y3 + H1
      IF (Y .GE. A) GO TO 81
      Y = A
      GO TO 90
   81 IF (Y .LE. B) GO TO 82
      Y = B
      GO TO 90
   82 CALL CERR(Y, F(Y), G(Y), PHI(Y), DEL, IERR, L, LP1, M, NP1, D)
      IF (IERR .NE. 0) RETURN
      Z = U*DEL
      IF (Z .LE. Z3) GO TO 83
      Y2 = Y3
      Y3 = Y
      Z2 = Z3
      Z3 = Z
      GO TO 80
   83 Y = (Z - Z3) + (Z2 - Z3)
      IF (Y .NE. ZERO) GO TO 84
      Y = Y3
      GO TO 90
   84 Y = HALF*(Y2 + Y3) + H1*(Z2 - Z3)/Y
C
   90 X(I) = Y
      CALL CERR(Y, F(Y), G(Y), PHI(Y), DEL, IERR, L, LP1, M, NP1, D)
      IF (IERR .NE. 0) RETURN
      ERR(I) = DEL
      U = -U
      IF (I .EQ. 1) GO TO 91
      IF (X(I) .LE. X(I-1)) GO TO 230
   91 Z = DABS(ERR(I))
      ERROR = DMAX1(ERROR, Z)
      IF (Z .GE. TEN) GO TO 240
      Y = DABS(XLB)
      ZZ = ONE
      IF (Y .NE. ZERO) ZZ = DABS(Z - Y)/Y
      IF (Z1 .LT. ZZ) Z1 = ZZ
   92 CONTINUE
C
C     SEARCH FOR AN EXTRA EXTREMAL POINT BETWEEN THE ENDPOINTS
C     OF THE INTERVAL AND THE CRITICAL POINTS
C
      IF (X(1) .LE. A) GO TO 110
      H1 = C0*(X(1) - A)
      U = ONE
      IF (XLB .GE. ZERO) U = -U
      Z3 = ZERO
      Y = A
      DO 100 I = 1,16
      CALL CERR(Y, F(Y), G(Y), PHI(Y), DEL, IERR, L, LP1, M, NP1, D)
      IF (IERR .NE. 0) RETURN
      Z = U*DEL
      IF (Z .LE. Z3) GO TO 100
      Z3 = Z
      Z2 = Y
  100 Y = Y + H1
      ERROR = DMAX1(ERROR, Z3)
      Z = DABS(XLB)
      IF (Z3 .LE. Z) GO TO 110
      I = NP1
      DO 101 II = 2,NP1
      ERR(I) = ERR(I-1)
      X(I) = X(I-1)
  101 I = I - 1
      X(1) = Z2
      ERR(1) = U*Z3
      GO TO 113
C
  110 IF (X(NP1) .GE. B) GO TO 120
      H1 = C0*(B - X(NP1))
      U = ONE
      IF (ERR(NP1) .GE. ZERO) U = -U
      Z3 = ZERO
      Y = B
      DO 111 I = 1,16
      CALL CERR(Y, F(Y), G(Y), PHI(Y), DEL, IERR, L, LP1, M, NP1, D)
      IF (IERR .NE. 0) RETURN
      Z = U*DEL
      IF (Z .LE. Z3) GO TO 111
      Z3 = Z
      Z2 = Y
  111 Y = Y - H1
      ERROR = DMAX1(ERROR, Z3)
      Z = DABS(XLB)
      IF (Z3 .LE. Z) GO TO 120
      DO 112 I = 1,N
      ERR(I) = ERR(I+1)
  112 X(I) = X(I+1)
      X(NP1) = Z2
      ERR(NP1) = U*Z3
  113 XLB = -XLB
      ZZ = ONE
      IF (Z .NE. ZERO) ZZ = DABS(Z3 - Z)/Z
      IF (Z1 .LT. ZZ) Z1 = ZZ
C
C     CHECK FOR CONVERGENCE
C
  120 IF (Z1 .LE. EPS) RETURN
C
C     SET UP FOR THE NEXT ITERATION
C
      IF (ITNO .GE. MXITER) GO TO 210
      SUM = ZERO
      SIGN = ONE
      DO 130 I = 1,NP1
      SUM = SUM + SIGN*ERR(I)
  130 SIGN = -SIGN
      XLB = SUM/DNP1
      ITNO = ITNO + 1
      GO TO 30
C     ------------------------------------------------------------------
C     ERROR RETURN
C     ------------------------------------------------------------------
C     INPUT ERROR
C
  200 IERR = 1
      RETURN
C
C     MXITER ITERATIONS WERE PERFORMED - MORE ITERATIONS ARE NEEDED
C
  210 IERR = 2
      RETURN
C
C     THE LINEAR EQUATIONS CANNOT BE SOLVED
C
  220 IF (ITNO .EQ. 1) GO TO 250
      IERR = 3
      RETURN
C
C     THE SEQUENCE OF CRITICAL POINTS IS NOT MONOTONICALLY INCREASING
C
  230 IERR = 4
      IF (I .LE. N) ERROR = OLDERR
      RETURN
C
C     IT APPEARS THAT THE ALGORITHM HAS FAILED TO CONVERGE
C     THERE MAY BE POLES IN THE RATIONAL APPROXIMATION
C
  240 IERR = 5
      RETURN
C
C     THE ROUTINE HAS COMPLETELY FAILED - THE RESULTS SHOULD BE IGNORED
C
  250 IERR = 6
      RETURN
      END
