      FUNCTION PR(X0, X, F, N, ERRMSG)

C     Polynomial approximant at X0
C     Limited to N = 10

      PARAMETER ( MAXN = 10 )
      REAL X0, X(N), F(N), FT(MAXN), H(MAXN)
      CHARACTER*80 ERRMSG

      ERRMSG = ' '

C     ------ Initialize arrays
      DO 1000 K = 1, N
         H(  K ) = X( K ) - X0
         FT( K ) = F( K )
 1000 CONTINUE

C     ------ Recursion for solution
      IF ( N .GE. 2) THEN
         DO 3000 I = 1, N - 1
            DO 2000 J = 1, N - I
               FT(J) = ( H( J + I ) * FT( J ) - H( J ) * FT( J + 1 ) ) /
     &                           ( H( J + I ) - H( J ) )
 2000       CONTINUE
 3000    CONTINUE
      ENDIF

      PR = FT( 1 )

      RETURN
      END
