      COMPLEX FUNCTION PADE( X0, X, F, N, ERRMSG )

C     COMPUTES THE [MU, NU] PADE APPROXIMANT AT X0, WHERE
C        MU = [M/2], NU = M - MU AND M = N - 1

C     FORMULAS TAKEN FROM BURLISCH AND STOER
C        NUMERISCHE MATHEMATIK 8:1-13(1966)

C     LIMITED TO N = 10 (OR SO)

C     N IS THE NUMBER OF POINTS USED.
C     X CONTAINS THE INDEPENDENT VARIABLE
C     F CONTAINS THE FUNCTIONAL VALUES AT EACH X
C     X0 IS THE POINT AT WHICH THE APPROXIMATION IS TO BE CALCULATED

C     NOTE MUCH WASTED SPACE IN THIS IMPLEMENTATION

C     MICHAEL PORTER, 11 SEPT 1985

      PARAMETER ( MAXN = 10, PREC = 1.0E-7 )
      COMPLEX X0, X( N ), F( N ), RDT, WDEN, CNUM, CDEN, 
     &    DT( MAXN, MAXN ), C( MAXN, MAXN ), H( MAXN ), XINT
      CHARACTER*80 ERRMSG

      ERRMSG = ' '

C     ------ Initialize arrays
      NEARST = 1
      HSMALL = 1.0E9

      DO K = 1, N
         H( K ) = X( K ) - X0

C        ------ Check for nearest abscissa
         IF ( ABS( H( K ) ) .LT. HSMALL ) THEN
            HSMALL = ABS( H( K ) )
            NEARST = K

C           ------ Compute x-interval size
            IF ( K .EQ. 1 ) THEN
               XINT = X( 2 ) - X( 1 )
            ELSE IF ( K .EQ. N ) THEN
               XINT = X( N ) - X( N - 1 )
            ELSE
               XINT = ( X( K + 1 ) - X( K - 1 ) ) / 2.0
            ENDIF

C           ------ Quick return if x lies "on" a data point
            IF ( HSMALL .LT. PREC * ABS( XINT ) ) THEN
               PADE = F( K )
               RETURN
            ENDIF
         ENDIF
      END DO

C     *** Recursion for solution ***

      DO M = 1, N
         DT( 1, M ) = F( M )
         C(  1, M ) = F( M )

         IF ( M .GE. 2 ) THEN

            DO K = 1, M - 1
               RDT  = H( M - K ) / H( M ) * DT( K, M - K )
               CNUM = C( K, M - K + 1 ) - DT( K, M - K )
               CDEN = RDT - C( K, M - K + 1 )

               IF ( ABS( CNUM ) .GE. 1.0E10 * ABS( CDEN ) ) THEN
C                  ERRMSG = 'ERROR: NEARLY-SINGULAR MATRIX IN PADE'
                  PADE = F( NEARST )
                  RETURN
               ENDIF

               WDEN = CNUM / CDEN
               DT( K + 1, M - K ) = C( K, M - K + 1 ) * WDEN
               C(  K + 1, M - K ) = RDT * WDEN
            END DO
         ENDIF
      END DO

C     ------ Sum side elts to get interpolate
      PADE = 0.0
      DO K = 1, N
         PADE = PADE + DT( K, N - K + 1 )
      END DO

      RETURN
      END
