function [divdif, p] = interpolation(x,y1)
% Die Funktion interpolation soll die Höchstkoeffizienten (= dividierte Differenzen)
% des Interpolationspolynom berechnen
% -------------------------------------------------------------------------
% EINGABE:
% x - Knoten, hier: x = [0 0.5 1 1.5 2]
% y - Datenpunkte, hier: y1 = [1 0 -1 0 1]
% AUSGABE:
% div - dividierte Differenzen
% p - Funktionswert an der Stelle y, hier: y2 = [0.25 0.75 1.25 1.75]
% -------------------------------------------------------------------------
x = [0 0.5 1 1.5 2]; 
y1 = [1 0 -1 0 1];

% Wir berechnen nun die dividierten Differenzen 

n = length(x);                      % gibt die Länge von x an, also die Anzahl der Knoten
divdif = zeros(1,n+1);              % erzeugt einen 1x(n+1) - Nullmatrix

for i = 1 : n+1                     % entspricht i = 0,...,n aus S. 19 vom Skript
    t = divdif(1);                  % setzt t gleich dem 1. Eintrag in der Nullmatrix
    divdif(1) = f(i);               % setzt diesen 1. Eintrag in der Nullmatrix gleich dem i-ten Funktionswert
    for k = 1 : i-1                 % entspricht k = 1,...,n aus S. 19 vom Skript
        Divdif = t;                 % setzt die divergierte Differenz mit einem xi mehr auf t
        t = divdif(k+1);            % setzt t gleich dem (k+1)ten Eintrag in der Nullmatrix
        divdif(k+1) = (divdif(k)-Divdif)./(x(i)-x(i-k));    % Formel für die dividierten Differenzen
    end
    d(i) = divdif(i);               % setzt d(i) gleich der i-ten divergierten Differenz
    disp([x(i) divdif(1:i)])        % zeigt die divergierten Differenzen an
end

% Nun wollen wir das Interpolationspolynom an der Stelle y ausrechnen. Wir
% benutzen dazu das Horner - Schema.

p = zeros(length(y));                   % 4x4 - Nullmatrix, da y hier 4 Einträge hat
for i = n:-1:1                        % entspricht i = 5,4,3,2,1
    p = (y - x(i)).*p + divdif(i);      % vgl. Formel S. 20 aus dem Skript
end
end