> restart;
> with(LinearAlgebra):
Newtoninterpol:=proc(x,f,z)
# x ist der Vektor der Stuetzstellen
# f ist der Vektor der Stuetzwerte
# z ist der Vektor der xKoordinaten, an denen der WErt des Interpolationspolynoms gesucht ist
local Newtonkoeff, y,p,m,n,j,k;
Newtonkoeff:=proc(x,f)
# x ist der Vektor der Stuetzszellen
# f ist der Vektor der Stuetzwerte
# option trace;
local y,n,i,k;
n:=Dimension(x);
y:=f;
for k from 1 to n-1 do
    for i from n by -1 to k+1 do
       y[i]:=evalf((y[i]-y[i-1])/(x[i]-x[i-k]));
    end do;
end do;
evalm(y);
end proc;
# Berechnung der Werte des Interplationspolynoms
y:=Newtonkoeff(x,f);
n:=Dimension(x);
m:=Dimension(z);
p:=Vector(m);
for j from 1 to m do
   p[j]:= y[n];
   for k from n-1 by (-1) to 1 do
       p[j]:=evalf(y[k]+(z[j]-x[k])*p[j]);
   end do;
end do;
evalm(p);
end proc:   
> n:=100:
  x:=Vector([0,1,3]);
> y:=Vector([1,3,2]);
  z:=Vector(n):
  a:=0; b:=3;
  for j from 1 to n do z[j]:=evalf(a+(j-1)*(b-a)/(n-1)); end do:
  p:=Newtoninterpol(x,y,z):
  with(plots):
  listplot([seq([z[t],p[t]],t=1..n)]);