	PROGRAM polint
	USE types

        integer, parameter :: ngx=15  ! size of the grid
        integer, parameter :: n=3     ! order of interpolation
	
        REAL(SP), DIMENSION(ngx) :: grid,coord 
        REAL(SP), DIMENSION(n) :: xa,ya
	REAL(SP) x
	REAL(SP) :: y,dy
	INTEGER(I4B) :: m,ns,i
	REAL(SP), DIMENSION(n) :: c,d,den,ho

!       initialize coordinates and function values 
        forall(i=1:ngx) coord(i)=i-1
        forall(i=1:ngx:2) grid(i)=3*coord(i)+1 !initalize every second
!                                          the rest is to be inerpolated


do i=1,ngx-n,2

        x=coord(i+1)
        xa(1:n)=coord(i:i+2*n-1:2)
        ya(1:n)=grid(i:i+2*n-1:2)

	c=ya
	d=ya
	ho=xa-x
	ns=1
	y=ya(ns)
	ns=ns-1
	do m=1,n-1
		den(1:n-m)=ho(1:n-m)-ho(1+m:n)
		den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
		d(1:n-m)=ho(1+m:n)*den(1:n-m)
		c(1:n-m)=ho(1:n-m)*den(1:n-m)
	        dy=c(ns+1)
		y=y+dy
	end do

write(6,*) 'interpolation polynomial'
write(6,*) xa
write(6,*) ya
write(6,*) 'x-value',x
write(6,*) 'solution',y
write(6,*) 'error',dy

enddo

	END program polint
