> restart;
> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> Newtonkoeff:=proc(x,f)
> #x ist der Vektor der Stuetzstellen
> #f ist der Vektor der Stuetzwerte
> local y,n,i,k;
> n:=vectdim(x);
> y:=vector(n);
> for i from 1 to n do y[i]:=f[i] od;
> 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]))
>    od
> od;
> evalm(y)
> end:
> 
> x:=vector(3,[0,1,3]):
> y:=vector(3,[1,3,2]):
> w:=Newtonkoeff(x,y);

                      w := [1, 2., -.8333333333]

> Newtoninterpol:=proc(x,f,z)
> # f ist der Vektor der Stuetzwerte
> # x ist der Vektor der Stuetzstellen
> # z ist der Vektor der Stuetzstellen, fuer den ein Interpolationswert gesucht ist
> local Newtonkoeff,y,p,m,n,j,k;
> Newtonkoeff:=proc(x,f)
> #x ist der Vektor der Stuetzstellen
> #f ist der Vektor der Stuetzwerte
> local y,n,i,k;
> n:=vectdim(x);
> y:=vector(n);
> for i from 1 to n do y[i]:=f[i] od;
> 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]))
>    od
> od;
> evalm(y)
> end;
> # Berechnung der Funktionswerte des Interpolationspolynoms
> y:=Newtonkoeff(x,f);
> n:=vectdim(y);
> m:=vectdim(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]) od
> od;
> evalm(p)
> end:
> 
> n:=100: z:=vector(n):
> a:=0: b:=3:
> for j from 1 to n do z[j]:=evalf(a+(j-1)*(b-a)/(n-1)) od:
> p:=Newtoninterpol(x,y,z):
> with(plots):
> listplot([seq([z[t],p[t]],t=1..n)]);

>