> restart;
with(LinearAlgebra):
Digits:=20;
#Neville Algorithmus
neville := proc(x,f,wert)
# x Vektor der Stuetzstellen
# f Vektor der Stuetzwerte
# wert ist die x-Koordinate an der das Interpolationspolynom ausgewertet werden soll
# option trace;
local y,n,j,k;
n:=Dimension(f);
y:=Vector(n);
y:=f;
for k from 1 to n-1 do
    for j from n by (-1) to k+1 do
      y[j]:= evalf(((wert- x[j-k])*y[j] + (x[j]-wert)*y[j-1])/(x[j]-x[j-k]));
    end do;
end do;
RETURN(y[n]);
end proc:

> 
> 
                               20
> x:=Vector([1.4,1.5,1.6,1.7]);
y:=Vector([0.9523,0.9661,0.9763,0.9838]);
w:=neville(x,y,1.43);
                 Vector[column](%id = 6399744)
                 Vector[column](%id = 6399864)
                     0.95687155000000000000
>