program INTERPOLATION
   implicit none
   integer, parameter , unit(n) :: N = 5
   real :: F,DF
   real, unit(x) :: X
   real, dimension (N) :: XI,FI
   data /XI/0.0,0.5,1.0,1.5,2.0
/FI/1.0,0.938470,0.765198,0.511828,0.223891
X = 0.9
call AITKEN(N,XI,FI,X,F,DF)
write (number = 6,'"(3F16.8)"') X,F,DF
end program INTERPOLATION
subroutine AITKEN(N,XI,FI,X,F,DF)
   integer, parameter , unit(n) :: NMAX = 21
   integer, intent(in)  :: N
   integer :: I,J
   real, intent(in)  :: X
   real, intent(out)  :: F,DF
   real :: X1,X2,F1,F2
   real, intent(in) , dimension (N) :: XI,FI
   real, dimension (NMAX) :: FT
if (N>NMAX) then
stop 'Dimension of the data is too large.'
   end if
do I = 1, N, 1
FT(I) = FI(I)
   end do
do I = 1, N-1, 1
do J = 1, N-I, 1
X1 = XI(J)
X2 = XI(J+I)
F1 = FT(J)
F2 = FT(J+1)
FT(J) = ((X-X1)/(X2-X1))*F2+((X-X2)/(X1-X2))*F1
      end do
   end do
F = FT(1)
DF = (ABS(F-F1)+ABS(F-F2))/2.0
end subroutine AITKEN