***NEWTON INTERPOLATION***************************************************
*     PROGRAM FOR NEWTON'S DIVIDED DFFERENCE INTERPOLATING POLYNOMIAL    *
**************************************************************************
*
      program NEWTON_INTERPOL
      implicit doubleprecision(a-h,o-z)
      parameter (md = 20, nd = 40)
      doubleprecision X(md),FX(nd,nd)
*
      print*, 'Enter number of data points, n:'
      read(*,*) n
      print*, 'Enter "X" data:'
      read(*,*) (X(i),i = 1,n)
      print*, 'Enter "FX" data:'
      read(*,*) (FX(i,1),i = 1,n)
      write(*,*)'------------------------------------------------------'
c
c     compute the finite divided differences
      m = n - 1
      do j = 1, n
         k = j + 1
         np = n - j
         do i = 1, np
            FX(i,k) = (FX(i+1,j) - FX(i,j))/(X(i+j)-X(i))
         enddo 
      enddo 
*
      write(6,*) 'The coefficients of interpolating polynomial,'
      write(6,*) ' F(x) = b0 + b1(x-x1) + b2(x-x1)(x-x2) + .... +bn(x-x1
     $)(x-x2)...(x-xn)   are:'
      write(6,*)
      write(6,*) (FX(1,j), j = 1,n)
*
      icontinue = 1
10    if (icontinue == 1) then
         write(6,*)
         write(6,*) 'Enter independent variable, where interpolation is 
     $desired, x:'
         read(*,*) xv
         fa = 1
         y = 0
         write(6,*)'----------------------------------------------------
     $---'
         write(6,101) 'PREDICTED VALUE','ERROR'
         do j = 1, n
            y = y + FX(1,j) *fa
            fa = fa *(xv - X(j))
            if(j >= n) go to 50
            jp = j + 1
            error = fa *FX(1,jp)
            write(6,102) j-1,'- Order',y,error
         enddo 
50       write(6,103) j-1,'- Order',y 
         write(6,*)'----------------------------------------------------
     $---'
         write(6,*) ' More Interpolation?  < 0 for NO / 1 for YES >:'
         read(*,*) icontinue
         go to 10
      endif
*
101   format('',24x,a15,8x,a5)
102   format('',3x,i1,a8,13x,f10.4,12x,f8.4) 
103   format('',3x,i1,a8,13x,f10.4)
*
      stop
      end