function [x,s]=myspline(a,b,f,n,m)

% m-file zur Uebung vom 29.5.09

% erstellt einen natuerlichen interpolierenden kubischen Spline s zur Funktion f 
% auf dem Intervall [a,b] mit n+1 aequidistanten Stuetzstellen, und wertet die 
% Splinefunktion an n*m+1 aequidistanten Stellen im Intervall [a,b] aus

%input: a,b reelle Zahlen; Intervallgrenzen 
%	f: String; zu interpolierende Funktion
% 	n,m: natuerliche Zahlen; n+1 ist die Anzahl der 
%	     aquidistanten Stuetzstellen xs_0...xs_n im Intervall [a,b]
%	     m+1 ist die Anzahl der auszuwertenden (aequidistanten) Stellen in jedem 
%	     Intervall [xs_i,xs_i+1]

%output: Vektor x der Laenge n*m+1 der auszuwertenden Stellen
%	 Vektor s der Laenge n*m+1 mit den Werten s(x)



try
% Erzeugen von n+1 aequidistanten Stuetzstellen
x_stuetz=linspace(a,b,n+1)';
y_stuetz=feval(f,x_stuetz);

% Vektor der aequidistanten auszuwertenden Stellen, je m im Intervall[x_i,x_i+1); Endpunkt nicht vergessen
x=linspace(a,b,m*n+1);
y=feval(f,x);


% Berechnung der Momente
h=x_stuetz(2:n+1)-x_stuetz(1:n);
M=momente(x_stuetz,y_stuetz);

% Berechnung der Koeffizienten
alpha=y_stuetz(1:n);
beta=(y_stuetz(2:n+1)-y_stuetz(1:n))./h-h/6.*(2*M(1:n)+M(2:n+1));
gamma=M(1:n)/2;
delta=1./(6*h).*(M(2:n+1)-M(1:n));


%s(x) berechnen
s=[];
for i=1:n
    % Spline erweitern um Spline auf [x_i,x_(i+1)]
    s=[s,alpha(i)+beta(i)*(x((i-1)*m+1:i*m)-x_stuetz(i))+gamma(i)*(x((i-1)*m+1:i*m)-x_stuetz(i)).^2+delta(i)*(x((i-1)*m+1:i*m)-x_stuetz(i)).^3];
end
s=[s,feval(f,b)]; % Endpunkt nicht vergessen

plot(x,s,x,y);
title('Runge Funktion und natuerlicher interpolierender Spline')
xlabel('x');
ylabel('y');
legend('Spline', 'Runge-Funktion')

catch
disp('Bitte Eingaben ueberpruefen')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function m=momente(x,y)
% berechnent die Momente m_0 .. m_n fuer gegebene Stuetzstellen und Funktionswerte für natuerliche RB
% input:    x:= Vektor der Länge n+1 mit den Stützstellen x_0...x_n
%           y:= Vektor der Länge n+1 mit den dazugehörigen Funktionswerten y_0...y_n
%           
% output:    m:= Vektor der Länge n+1 mit den Momenten m_0 ... m_n bei natuerlichen Randbedingungen


n=length(x);
h=x(2:n)-x(1:n-1);
deltay=y(2:n)-y(1:n-1);

% Erstelle S und c zur Loesung des Gleichungssystem aus der VL zur Bestimmung von m_1...m_n-1
c=6*deltay(2:n-1)./h(2:n-1)-6*deltay(1:n-2)./h(1:n-2);
S=diag(2*(h(1:n-2)+h(2:n-1)))+diag(h(2:n-2),1)+diag(h(2:n-2),-1);
m=S\c;

% natuerliche Randbedingungen beruecksichtigen (m_0=m_n=0)
m=[0;m;0];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=runge25(x)
% unsere Beispielfunktion
y=1./(1+25*x.^2);