% Anzahl der Daten
NumberOfData= 35;
% Stützstellen zum Interpolieren
SamplingPoints = linspace(-1,1,NumberOfData);
%Runge-Funktion
f = @(t) 1./(1+25.*t.*t);
% Funktionswerte berechnen
y = f(SamplingPoints);

% Punkte an denen später interpoliert werden soll
TestPoints = SamplingPoints;

DivDif = zeros(NumberOfData,NumberOfData);
DivDif(:,1) = y(:);
for j=2:NumberOfData
    for i=1:NumberOfData-j+1
        DivDif(i,j) = (DivDif(i+1,j-1)-DivDif(i,j-1))/(SamplingPoints(i+j-1)-SamplingPoints(i));
    end
end

xt = ones(1,NumberOfData);
yInterpol = DivDif(1,1);
for j=1:NumberOfData-1
    xt = xt.*(TestPoints-SamplingPoints(j));
    yInterpol = yInterpol+DivDif(1,j+1).*xt;
end


% links
plot(SamplingPoints,y,'r');
hold on
plot(SamplingPoints,yInterpol,'.k');