function [a] = newton_interpolation(x,f)

% newton_interpolation berechnet das Newton-Interpolationspolynom
% mit Hilfe von dividierten Differenzen.
% Aufruf:
%   [a] = newton_interpolation(x,f)
% x : St"utzstellen
% f : St"utzwerte
% a : Ausgabe, die Koeffizienten des Polynoms
% Die Matrix der dividierten Differenzen wird im Feld d gespeichert.

% berechne Dimension der Vektoren
[n,m] = size(x);
[nf,mf] = size(f);

% konvertiere, falls n"otig, die Vektoren zu Spaltenvektoren
if n==1
  x = x';
  n = m;
end

if nf==1
  f = f';
  nf = mf;
end

% kontrolliere, ob gleiche Anzahl von Werten vorhanden ist
if n ~= nf
  error('Die Dimensionen von x und f passen nicht')
end

% erniedrige die Dimension, damit die Formeln die gleiche Gestalt
% wie in der Vorlesung haben (Vorlesung: Indizes von 0 bis n;
% hier, Indizes von 1 bis n)
n = n-1;

% erzeuge die Matrix f"ur die dividierten Differenzen
d = zeros(n+1,n+1);

% belege erste Spalte mit der rechten Seite
d(:,1) = f;

% berechne die dividierten Differnzen
for k=2:n+1
  for i=k:n+1
    d(i,k) = (d(i,k-1)-d(i-1,k-1))/(x(i)-x(i-k+1));
  end
end

% belege den Vektor a mit der Hauptdiagonalen von d
a = diag(d);

% zeichne Bild des Interpolationspolynoms
% berechne Werte f"ur das Interpolationspolynom
steps = 1000;
h = (x(n+1)-x(1))/steps;

for i=1:steps+1;
  % Abszisse
  xp(i) = x(1) + (i-1)*h;
  % Verwende Formel (7) aus der Vorlesung
  % Initialisiere Funktionswert und Knotenpolynom
  fp(i) = 0;
  omega = 1;
  % berechne Funktionswert
  for j = 1:n+1
    fp(i) = fp(i) + a(j)*omega;
    omega = omega * (xp(i)-x(j));
  end
end

clf

plot(xp,fp,'-b',x,f,':r');





