      SUBROUTINE CHEBY ( NF, NPL, NPLMAX, N2, FUNCTN, X, FXJ, GC )

c*********************************************************************72
c
C  SIMULTANEOUS CHEBYSHEV ANALYSIS OF NF FUNCTIONS
C  COMPUTES A MATRIX, X, CONTAINING ONE CHEBYSHEV SERIES PER
C  COLUMN FOR A GIVEN NUMBER OF FUNCTIONS, NF.  INPUT NFL,
C  THE NUMBER OF TERMS IN ALL SERIES, NPLMAX, THE ROW
C  DIMENSION OF X IN THE CALLING PROGRAM (MUST BE.GE.NPL),
C  N2, DIMENSION OF GC (MUST BE.GE.2*(NPL-1)), AND FUNCTN,
C  THE NAME OF USER SUBROUTINE WHICH DEFINES THE NF
C  FUNCTIONS.  FXJ AND GC ARE WORK SPACE.
C  AN EXAMPLE OF SUCH A SUBROUTINE IS AS FOLLOWS
C  SUBROUTINE FUNCTN ( A, VAL )
C  DOUBLE PRECISION A, VAL(2)
C  VAL(1) = DSIN(A)
C  VAL(2) = DCOS(A)
C  RETURN
C  END
c
      DOUBLE PRECISION X(NPLMAX,NF), FXJ(NF), GC(N2), XJ,
     &  FK, PEN, FAC

      DO K = 1, NPL
        DO J = 1, NF
          X(K,J) = 0.D0
        end do
      end do
      N = NPL - 1
      ENN = N
      PEN = 3.1415926535897932D0 / ENN
      DO 30 K = 1, N2
        FK = K - 1
        GC(K) = DCOS ( FK * PEN )
30    CONTINUE
      DO 80 J = 1, NPL
        XJ = GC(J)
        CALL FUNCTN ( XJ, FXJ )
        IF ( J .NE. 1 .AND. J .NE. NPL ) GO TO 50
        DO 40 K = 1, NF
          FXJ(K) = .5D0 * FXJ(K)
40      CONTINUE
50      DO 70 L = 1, NPL
          LM = MOD ( ( L - 1 ) * ( J - 1 ), N2 ) + 1
          DO 60 K = 1, NF
            X(L,K) = X(L,K) + FXJ(K) * GC(LM)
60        CONTINUE
70      CONTINUE
80    CONTINUE
      FAC = 2.0D0 / ENN
      DO 100 K = 1, NPL
        DO 90 J = 1, NF
          X(K,J) = FAC * X(K,J)
90      CONTINUE
100   CONTINUE
      RETURN
      END
      SUBROUTINE MLTPLY ( XX, X2, NPL, X3 )

c*********************************************************************72
c
C  MULTIPLIES TWO GIVEN CHEBYSHEV SERIES, XX AND X2, WITH
C  NPL TERMS TO PRODUCE AN OUTPUT CHEBYSHEV SERIES, X3.
c
      DOUBLE PRECISION XX(NPL), X2(NPL), X3(NPL), EX

      DO 10 K = 1, NPL
        X3(K) = 0.0D0
10    CONTINUE
      DO 30 K = 1, NPL
        EX = 0.0D0
        MM = NPL - K + 1
        DO 20 M = 1, MM
          L = M + K - 1
          EX = EX + XX(M) * X2(L) + XX(L) * X2(M)
20      CONTINUE
        X3(K) = 0.5D0 * EX
30    CONTINUE
      X3(1) = X3(1) - 0.5D0 * XX(1) * X2(1)
      DO 50 K = 3, NPL
        EX = 0.0D0
        MM = K - 1
        DO 40 M = 2, MM
          L = K - M + 1
          EX = EX + XX(M) * X2(L)
40      CONTINUE
        X3(K) = 0.5D0 * EX + X3(K)
50    CONTINUE
      RETURN
      END
      SUBROUTINE ECHEB ( X, COEF, NPL, FX )

c*********************************************************************72
c
C  EVALUATES THE VALUE FX(X) OF A GIVEN CHEBYSHEV SERIES,
C  COEF, WITH NPL TERMS AT A GIVEN VALUE OF X BETWEEN
C  -1. AND 1.
c
      DOUBLE PRECISION COEF(NPL), X, FX, BR, BRPP, BRP2

      BR = 0.0D0
      BRPP = 0.0D0
      DO 10 K = 1, NPL
        J = NPL - K + 1
        BRP2 = BRPP
        BRPP = BR
        BR = 2.0D00 * X * BRPP - BRP2 + COEF(J)
10    CONTINUE
      FX = 0.5D0 * ( BR - BRP2 )
      RETURN
      END
      SUBROUTINE EDCHEB ( X, COEF, NPL, FX )

c*********************************************************************72
c
C  EVALUATES THE VALUE FX(X) OF THE DERIVATIVE OF A
C  CHEBYSHEV SERIES, COEF, WITH NPL TERMS AT A GIVEN
C  VALUE OF X BETWEEN -1. AND 1.
c
      DOUBLE PRECISION COEF(NPL), X, FX, XJP2, XJPL, XJ, BJP2,
     &  BJPL, BJ, BF, DJ

      XJP2 = 0.0D0
      XJPL = 0.0D0
      BJP2 = 0.0D0
      BJPL = 0.0D0
      N = NPL - 1
      DO 10 K = 1, N
        J = NPL - K
        DJ = J
        XJ = 2.D0 * COEF(J+1) * DJ + XJP2
        BJ = 2.D0 * X * BJPL - BJP2 + XJ
        BF = BJP2
        BJP2 = BJPL
        BJPL = BJ
        XJP2 = XJPL
        XJPL = XJ
10    CONTINUE
      FX = .5D0 * ( BJ - BF )
      RETURN
      END
      SUBROUTINE DFRNT ( XX, NPL, X2 )

c*********************************************************************72
c
C  DFRNT COMPUTES THE DERIVATIVE CHEBYSHEV SERIES, X2, OF A GIVEN
C  CHEBYSHEV SERIES, XX, WITH NPL TERMS.
C  TO REPLACE A SERIES X BY ITS DERIVATIVE, USE
C  CALL DFRNT ( X, NPL, X )
c
      DOUBLE PRECISION XX(NPL), XXN, XXL, DL, X2(NPL)

      DN = NPL - 1
      XXN = XX(NPL-1)
      X2(NPL-1) = 2.D0 * XX(NPL) * DN
      X2(NPL) = 0.D0
      DO 10 K = 3, NPL
        L = NPL - K + 1
        DL = L
        XXL = XX(L)
        X2(L) = X2(L+2) + 2.D0 * XXN * DL
        XXN = XXL
10    CONTINUE
      RETURN
      END
      SUBROUTINE NTGRT ( XX, NPL, X2 )

c*********************************************************************72
c
C  COMPUTES THE INTEGRAL CHEBYSHEV SERIES, X2, OF A GIVEN
C  CHEBYSHEV SERIES, XX, WITH NPL TERMS.
C  TO REPLACE A SERIES X BY ITS INTEGRAL, USE
C  CALL NTGRT ( X, NPL, X )
c
      DOUBLE PRECISION XX(NPL), XPR, TERM, DK, X2(NPL)

      XPR = XX(1)
      X2(1) = 0.0D0
      N = NPL - 1
      DO 10 K = 2, N
        DK = K - 1
        TERM = ( XPR - XX(K+1) ) / ( 2.D0 * DK )
        XPR = XX(K)
        X2(K) = TERM
10    CONTINUE
      DK = N
      X2(NPL) = XPR / ( 2.D0 * DK )
      RETURN
      END
      SUBROUTINE INVERT ( X, XX, NPL, NET, XNVSE, WW, W2 )

c*********************************************************************72
c
C  COMPUTES THE INVERSE CHEBYSHEV SERIES, XNVSE, GIVEN A
C  CHEBYSHEV SERIES, X, A FIRST APPROXIMATION CHEBYSHEV
C  SERIES, XX, WITH NPL TERMS, AND THE NUMBER OF
C  ITERATIONS, NET.  THE SUBROUTINE USES THE EULER METHOD
C  AND COMPUTES ALL POWERS EPS**K UP TO K=2**(NET+1),
C  WHERE EPS = 1 - X * ( XX INVERSE ).  WW AND W2 ARE WORK SPACE.
C  SUBROUTINES USED - MLTPLY
c
      DOUBLE PRECISION X(NPL), XX(NPL), XNVSE(NPL), WW(NPL),
     &  W2(NPL)

      CALL MLTPLY ( X, XX, NPL, WW )
      WW(1) = 2.D0 - WW(1)
      DO 10 K = 2, NPL
        WW(K) = -WW(K)
10    CONTINUE
      CALL MLTPLY ( WW, WW, NPL, W2 )
      WW(1) = 2.D0 * WW(1)
      DO 40 K = 1, NET
        CALL MLTPLY ( WW, W2, NPL, XNVSE )
        DO 20 J = 1, NPL
          WW(J) = WW(J) + XNVSE(J)
20      CONTINUE
        CALL MLTPLY ( W2, W2, NPL, XNVSE )
        DO 30 J = 1, NPL
          W2(J) = XNVSE(J)
30      CONTINUE
40    CONTINUE
      CALL MLTPLY ( WW, XX, NPL, XNVSE )
      RETURN
      END
      SUBROUTINE BINOM ( X, XX, NPL, M, NT, XA, WW, W2, W3 )

c*********************************************************************72
c
C  COMPUTES THE BINOMIAL EXPANSION SERIES, XA, FOR (-1/M)
C  POWER OF A GIVEN CHEBYSHEV SERIES, X, WITH NPL TERMS,
C  WHERE M IS A POSITIVE INTEGER.  XX IS A GIVEN INITIAL
C  APPROXIMATION TO X**(-1/M).  NT IS A GIVEN NUMBER OF
C  TERMS IN BINOMIAL SERIES.  WW, W2 AND W3 ARE WORK SPACE.
C  SUBROUTINES USED - MLTPLY.
c
      DOUBLE PRECISION X(NPL), XX(NPL), XA(NPL), WW(NPL),
     &  W2(NPL), W3(NPL), ALFA, COEF, DM, DKMM, DKM2

      DM = M
      ALFA = -1.D0 / DM
      DO 10 J = 1, NPL
        WW(J) = X(J)
10    CONTINUE
      DO 30 K = 1, M
        CALL MLTPLY ( WW, XX, NPL, W2 )
        DO 20 J = 1, NPL
          WW(J) = W2(J)
20      CONTINUE
30    CONTINUE
      WW(1) = WW(1) - 2.D0
      XA(1) = 2.D0
      DO 40 J = 2, NPL
        XA(J) = 0.0D0
        W3(J) = 0.D0
40    CONTINUE
      W3(1) = 2.D0
      DO 60 K = 2, NT
        DKMM = K - 1
        DKM2 = K - 2
        COEF = ( ALFA - DKM2 ) / DKMM
        CALL MLTPLY ( W3, WW, NPL, W2 )
        DO 50 J = 1, NPL
          W3(J) = W2(J) * COEF
          XA(J) = XA(J) + W3(J)
50      CONTINUE
60    CONTINUE
      CALL MLTPLY ( XA, XX, NPL, W2 )
      DO 70 J = 1, NPL
        XA(J) = W2(J)
70    CONTINUE
      RETURN
      END
      SUBROUTINE XALFA2 ( X, XX, NPL, M, MAXET, EPSLN, NET, WW,
     &  W2 )

c*********************************************************************72
c
C  REPLACES A GIVEN INITIAL APPROXIMATION CHEBYSHEV SERIES,
C  XX, BY A GIVEN CHEBYSHEV SERIES, X, WITH NPL TERMS,
C  RAISED TO THE (-1/M) POWER, WHERE M IS AN INTEGER.
C  INPUT MAXET, MAXIMUM ALLOWED NUMBER OF ITERATIONS, AND
C  EPSLN, REQUIRED PRECISION EPSILON.  OUTPUT ARGUMENT,
C  NET, IS NUMBER OF ITERATIONS PERFORMED.  IF MAXET=NET,
C  MAY BE DIVERGENCE.  WW AND W2 ARE WORK SPACE.
C  CONVERGENCE IS QUADRATIC.
C  SUBROUTINES USED - MULTIPLY.
c
      DOUBLE PRECISION X(NPL), XX(NPL), WW(NPL), W2(NPL),
     & EPSLN, DALFA, DM, S, TDMM

      DM = M
      DALFA = 1.D0 / DM
      TDMM = 2.D0 * ( DM + 1.D0 )
      DO 60 JX = 1, MAXET
        DO 10 L = 1, NPL
          WW(L) = X(L)
10      CONTINUE
        DO 30 K = 1, M
          CALL MLTPLY ( WW, XX, NPL, W2 )
          DO 20 L = 1, NPL
            WW(L) = W2(L)
20        CONTINUE
30      CONTINUE
        S = -2.D0
        DO 40 L = 1, NPL
          S = S + DABS ( WW(L) )
          WW(L) = -WW(L)
40      CONTINUE
        WW(1) = WW(1) + TDMM
        CALL MLTPLY ( WW, XX, NPL, W2 )
        DO 50 L = 1, NPL
          XX(L) = W2(L) * DALFA
50      CONTINUE
        NET = JX
        IF ( DABS ( S ) .LT. EPSLN ) RETURN
60    CONTINUE
      RETURN
      END
      SUBROUTINE XALFA3 ( X, XX, NPL, M, MAXET, EPSLN, NET, WW,
     &  W2 )

c*********************************************************************72
c
C  REPLACES A GIVEN INITIAL APPROXIMATION CHEBYSHEV SERIES,
C  XX, BY A GIVEN CHEBYSHEV SERIES, X, WITH NPL TERMS,
C  RAISED TO THE (-1/M) POWER, WHERE M IS AN INTEGER.
C  INPUT MAXET, MAXIMUM ALLOWED NUMBER OF ITERATIONS, AND
C  EPSLN, REQUIRED PRECISION EPSILON.  OUTPUT ARGUMENT
C  NET, IS NUMBER OF ITERATIONS PERFORMED.  IF MAXET = NET,
C  MAY BE DIVERGENCE.  WW AND W2 ARE WORK SPACE.
C  CONVERGENCE IS OF ORDER THREE.
C  SUBROUTINES USED - MLTPLY
c
      DOUBLE PRECISION X(NPL), XX(NPL), WW(NPL), W2(NPL),
     &  EPSLN, DALFA, DM, S, TDMM, P5DML

      DM = M
      DALFA = 1.D0 / DM
      TDMM = 2.D0 * ( DM + 1.D0 )
      P5DML = .5D0 * ( DM + 1.D0 )
      DO 90 JX = 1, MAXET
        DO 10 L = 1, NPL
          WW(L) = X(L)
10      CONTINUE
        DO 30 K = 1, M
          CALL MLTPLY ( WW, XX, NPL, W2 )
          DO 20 L = 1, NPL
            WW(L) = W2(L)
20        CONTINUE
30      CONTINUE
        S = -2.D0
        DO 40 L = 1, NPL
          S = S + DABS ( WW(L) )
40      CONTINUE
        WW(1) = WW(1) - 2.D0
        DO 50 L = 1, NPL
          WW(L) = WW(L) * DALFA
50      CONTINUE
        CALL MLTPLY ( WW, WW, NPL, W2 )
        DO 60 L = 1, NPL
          WW(L) = -WW(L)
          W2(L) = W2(L) * P5DML
60      CONTINUE
        WW(1) = WW(1) + 2.D0
        DO 70 L = 1, NPL
          W2(L) = W2(L) + WW(L)
70      CONTINUE
        CALL MLTPLY ( W2, XX, NPL, WW )
        DO 80 L = 1, NPL
          XX(L) = WW(L)
80      CONTINUE
        NET = JX
        IF ( DABS ( S ) .LT. EPSLN ) RETURN
90    CONTINUE

      RETURN
      END
