Сначала я вычислил нижний треугольник и значения главной диагонали по формуле: f(i,j)=(f(i,j-1))-(f(i-1,j-1))/(x(i)- x(i-j), я сделал матрицу (6,6), чтобы я мог получить из нее основную диагональ, значения которой определяются с помощью b в коде и показывается как «Cosntante» в Promt.
Затем я выполнил интерполяцию, определенную как p в коде и показанную как "El valor aproximado de f () es:" в запросе.
program p4
use ISO_C_Binding
implicit none
real :: x(6), y(6), b(6), xi, p
integer :: i, j, n, m, k
real, dimension(6,6) :: f=c_float
f=0
x=[1.6, 2.0, 2.5, 3.2, 4.0, 4.5]
y=[2, 8, 14, 15, 8, 2]
n=SIZE(x)
do i=1,n
f(i,1)=y(i)
end do
do i=1,n
f(i,2:n)=0
end do
do j=2,n
do i=j,n
f(i,j)=(f(i,j-1))-(f(i-1,j-1))/(x(i)- x(i-j))
end do
end do
do i=1,n
write(*,*) F(i,:)
end do
do i=1,n
b(i)=f(i,i)
write(*,*) "Constante b[",i,"]=",b(i)
end do
xi=2.8
m=size(x)-1
p=b(m)
do k=2,m+1
p=(b(m-k))+(xi-((x(m-k))*p))
end do
write(*,*) "El valor aproximado de f(",xi,") es:",p
pause
end program p4
она показывает мне следующие значения (это значения, которые я получаю в ForceFortran): 
И это должно показать мне (это значения, которые я получаю в Python):
Что мне делать?