function MyLGS_Spline_N()
% Berechnet einen kubischen C^2-Spline durch die N Punkte P
% und liefert fuer N-1 Polynome p_j auf N-1 Intervallen
% die Koeffizienten a(*,J). Dann ist der Spline auf jedem 
% Intervall j definiert durch
%
%   p_j(x) = a(4,j)*x^3 + a(3,j)*x^2 + a(2,j)*x + a(1,j)
%
% Eingabe: nichts
% Ausgabe: ein plot
% Aufruf: >> MyLGS_Spline_N
%
close all

P = [[1 2];[2 -1];[3 1];[4 0];[6 1]];
N = length(P);

a = SplineKoef(P);

step = 1;
DIM  = 100;
for j=1:N-1
    h = (P(j+1,1)-P(j,1))/(DIM-1);
    for i=1:DIM
        x(step) = P(j,1) + (i-1)*h;
        y(step) = a(4,j)*x(step)^3+a(3,j)*x(step)^2+a(2,j)*x(step)+a(1,j);
        step = step + 1;
    end
end

plot(x,y,'b-',P(:,1),P(:,2),'mo')
grid on


end


function a=SplineKoef(P)

x = P(:,1); y = P(:,2);
N = length(P);
DIM = (N-1)*4;
A = zeros(DIM,DIM);

A(1,1:2) = [6*x(1) 2];                 
yy(1) = 0;
A(2,1:4) = [x(1)^3 x(1)^2 x(1) 1]; 
yy(2) = y(1);

for j=2:N-1
    A(4*(j-2)+3,4*(j-2)+1:4*(j-2)+4) = [x(j)^3 x(j)^2 x(j) 1];    yy(4*(j-2)+3) = y(j);
    A(4*(j-2)+4,4*(j-2)+5:4*(j-2)+8) = [x(j)^3 x(j)^2 x(j) 1];    yy(4*(j-2)+4) = y(j);
    A(4*(j-2)+5,4*(j-2)+1:4*(j-2)+4) = [3*x(j)^2 2*x(j) 1 0];     yy(4*(j-2)+5) = 0;
    A(4*(j-2)+6,4*(j-2)+5:4*(j-2)+8) = [-3*x(j)^2 -2*x(j) -1 0];  yy(4*(j-2)+6) = 0;
end

A(4*(N-1)-1,4*(N-1)-3:4*(N-1)) = [x(N)^3 x(N)^2 x(N) 1];
yy(4*(N-1)-1) = y(N);
A(4*(N-1),4*(N-1)-3:4*(N-1)) = [6*x(N) 2 0 0];
yy(4*(N-1)) = 0;

xx  = A\yy' %MyLR(A,yy)

for j=1:N-1
    a(4,j) = xx(4*(j-1)+1);
    a(3,j) = xx(4*(j-1)+2);
    a(2,j) = xx(4*(j-1)+3);
    a(1,j) = xx(4*(j-1)+4);
end

end

