	PROGRAM polint
	USE types

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

!       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

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

!forall(i=1:ngx-n:2) c(i,:)=ya(i,:)
forall(i=1:ngx-n:2) c(i,:)=grid(i:i+2*n-1:2)
write(6,*) c
forall(i=1:ngx-n:2) d(i,:)=grid(i:i+2*n-1:2)
forall(i=1:ngx-n:2) ho(i,:)=xa(i,:)-x(i)
forall(i=1:ngx-n:2) y(i)=ya(i,1)

forall(i=1:ngx-n:2) den(i,1)=ho(i,1)-ho(i,2)
forall(i=1:ngx-n:2) den(i,1)=(c(i,2)-d(i,1))/den(i,1)
!forall(i=1:ngx-n:2) d(i,1)=ho(i,2)*den(i,1)
forall(i=1:ngx-n:2) c(i,1)=ho(i,1)*den(i,1)
forall(i=1:ngx-n:2) y(i)=y(i)+c(i,1)

do i=1,ngx-n,2

write(6,*) 'interpolation polynomial'
write(6,*) xa(i,:)
write(6,*) grid(i:i+2*n-1:2)
write(6,*) 'x-value',x(i)
write(6,*) 'solution',y(i)

enddo

	END program polint
