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
 close all

% berechne Dimensionen f"ur das Bild
x_min = 1e10;
x_max = -1e10;
f_min = 1e10;
f_max = -1e10;

for i=1:n+1;
  if (x(i) < x_min)
    x_min = x(i);
  end
  if (x(i) > x_max)
    x_max = x(i);
  end
  if (f(i) < f_min)
    f_min = f(i);
  end
  if (f(i) > f_max)
    f_max = f(i);
  end
end

plot(x,f,'x');

axis([x_min-1,x_max+1,f_min-1,f_max+1])

hold on;

% berechne Werte f"ur das Interpolationspolynom
steps = 1000;
h = (x_max-x_min+2)/steps;

for i=1:steps+1;
  % Abszisse
  xp(i) = x_min-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

plot(xp,fp,'-r');





