function yi=interpNA(x,y,xi)
% Neville-Aitken Polynominterpolation
n=length(x)
ni=length(xi)
mati=zeros(ni,n);
for i=1:n
    mati(:,i)=y(i);   
end
% Invariante: mati(:,i) interpoliert in xi(i-k),...xi(i).
% Gilt zuerst fuer k=0, dann fuer alle k
for k=1:n-1
    for i=n:-1:k+1
        mati(:,i)=((xi-x(i-k)).*mati(:,i)+(x(i)-xi).*mati(:,i-1))...
            /(x(i)-x(i-k));   
    end 
end
% Ergebnis fuer k=n-1, i=n.
yi=mati(:,n);