%Interpolation nach Newton. Rungebeispiel mit aequidistanten Stuetzstellen.
clc
NumberOfData= 25;
SamplingPoints = linspace(-1,1,NumberOfData); %Generieren der n Stützstellen.
f = @(x) 1./(1+25.*x.*x); %Runge-Beispiel

TargetValues = f(SamplingPoints); %Dividierte Differenzen

for k=1:NumberOfData
    h(k) = TargetValues(k);
    if k>1
        for i=k-1:-1:1
            h(i) = (h(i+1)-h(i)) ./(SamplingPoints(k)-SamplingPoints(:,i));
        end
    end
    p(k) = h(1);
end

% zurück rechnen
n = NumberOfData;
vls = p(:,n)*ones(1,NumberOfData);
for i = NumberOfData-1:-1:1
    vls = vls.*(SamplingPoints-SamplingPoints(i)*ones(1,NumberOfData))+p(:,i)*ones(1,NumberOfData);
end


% links
plot(SamplingPoints,TargetValues,'r');
hold on
plot(SamplingPoints,vls,'.k');
legend('original Fkt','newton. Interp.');