function r = cubic_spline(x,f)

% cubic_spline berechnet den kubischen Spline durch die
% St"utzstellen (x(i),f(i)), i=1,...,n
% Aufruf r = cubic_spline(x,f)
% die Vorgehensweise aus der Vorlesung wird benutzt
% beachte: in Vorlesung n+1 St"utzstellen, hier nur n
 

% kontrolliere, ob die Dimensionen stimmen
[n,m] = size(x);
[n1,m1] = size(f);
if n ~= n1
   r = -1;
   error('Die Dimensionen der Felder stimmen nicht !');
end

% generiere System f"ur die Koeffizienten M_i
% belege Diagonale mit Zweien
A=2*eye(n,n);
for i=2:n-1
  % lambda_i
  A(i,i+1) = (x(i+1) - x(i))/(x(i+1)-x(i-1));
  % mu_i
  A(i,i-1) = (x(i)-x(i-1))/(x(i+1)-x(i-1));
  % rechte Seite d_i
  d(i) = 6*((f(i+1)-f(i))/(x(i+1)-x(i))-(f(i)-f(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1));
end
% die Werte f\"ur den Rand sind ohnehin auf Null gesetzt,
% f"ur die vorgegebene Bedingung also richtig, nur der letzte Wert der rechten
% Seite muss noch gesetzt werden, damit der Vektor die richtige Dimension besitzt
d(n) = 0;

M = A\d';

% berechne die Werte c_i und \tilde c_i

for i=2:n
  tilde_c(i-1,1) = f(i-1) - M(i-1)*(x(i)-x(i-1))^2/6;
  c(i-1,1) = (f(i) - f(i-1))/(x(i)-x(i-1)) - (x(i)-x(i-1))*(M(i)-M(i-1))/6;
end


% berechne Vektoren f"ur das Plotten
for i=2:n
  % in jedem Teilintervall 10 Funktionswerte
  h = (x(i) - x(i-1))/10;
  for j=1:11
    index = (i-2)*10+j;
    xp(index) = x(i-1) + (j-1)*h;
    % nutze Formel aus der Vorlesung, das Argument ist xp(index)	
    fp(index) = M(i-1) * (x(i)-xp(index))^3/(6*(x(i)-x(i-1)));
    fp(index) = fp(index) + M(i) * (xp(index)-x(i-1))^3/(6*(x(i)-x(i-1)));
    fp(index) = fp(index) + c(i-1)*(xp(index)-x(i-1))+tilde_c(i-1);
  end
end
clf
plot(xp,fp,'-b',x,f,':r');


