subroutine divdifftable(m, y, d)
!
!     Given a table of m+3 ordinates y of a function at the m+3
!     integer abscissae -1, 0, 1, ..., m+1, this subroutine
!     computes the Newton divided differences table d for the
!     interpolating 3rd-order polynomials on the abscissae
!     intervals [0..1], ..., [m-1..m].
!
!     B. Knapp, 2002-03-19
!
!     Reference: J. Vandergraft, Introduction to Numerical Computations,
!     2nd ed, Academic Press, 1983, algorithm 4.10 (p. 96).
!
      implicit none
!
!     Input
      integer*4 m
      real*8 y(-1:m+1)
!
!     Output
      real*8 d(4,0:m-1)
!
!     Local
      integer*4 i, j, k
!
!     Functions, subroutines
      intrinsic dble
!
      do j=0,m-1
         do i=1,4
            d(i,j) = y(j+i-2)
         enddo
         do k=1,3
            do i=4,k+1,-1
               d(i,j) = (d(i,j)-d(i-1,j))/dble(k)
            enddo
         enddo
      enddo
      return
      end