subroutine divdiffeval(m, d, x, y, dy)
!
!     Given a table of 3rd-order divided differences d (computed by
!     subroutine divdifftable), this routine evaluates the inter-
!     polating polynomial at the given abscissa x, and returns the 
!     interpolated value y and its truncation error dy.  The tabular 
!     abscissae used in the construction of the divided differences 
!     table must have been the integers -1, 0, 1, ..., m+1.
!
!     B. Knapp, 2002-03-19
!
!     Reference: J. Vandergraft, Introduction to Numerical Computations,
!     2nd ed., Academic Press 1983, Algorithm 4.12 (pp 98-103).
!
!     Input
      integer*4 m
      real*8 d(4,0:m-1), x
!
!     Output
      real*8 y, dy
!
!     Local
      integer*4 j, k
      real*8 z
!
!     Functions, subroutines
      intrinsic int, dble
!
      k = int(x)
      z = x-dble(k-1)
      y = d(1,k)+d(2,k)*z
      do j=2,3
         z = z*(x-dble(k+j-2))
         dy = d(j+1,k)*z
         y = y+dy
      enddo
      return
      end