function ex09()
xi = [1, 2, 3, 4, 5]';
yi = [1, 2, 1, 1, 4]';

%% Polynom
syms x;
P = (x-xi(2))*(x-xi(3))*(x-xi(4))*(x-xi(5))/...
    ((xi(1)-xi(2))*(xi(1)-xi(3))*(xi(1)-xi(4))*(xi(1)-xi(5))) * yi(1)...
+(x-xi(1))*(x-xi(3))*(x-xi(4))*(x-xi(5))/...
    ((xi(2)-xi(1))*(xi(2)-xi(3))*(xi(2)-xi(4))*(xi(2)-xi(5))) * yi(2)...
+(x-xi(1))*(x-xi(2))*(x-xi(4))*(x-xi(5))/...
    ((xi(3)-xi(1))*(xi(3)-xi(2))*(xi(3)-xi(4))*(xi(3)-xi(5))) * yi(3)...
+(x-xi(1))*(x-xi(2))*(x-xi(3))*(x-xi(5))/...
    ((xi(4)-xi(1))*(xi(4)-xi(2))*(xi(4)-xi(3))*(xi(4)-xi(5))) * yi(4)...
+(x-xi(1))*(x-xi(2))*(x-xi(3))*(x-xi(4))/...
    ((xi(5)-xi(1))*(xi(5)-xi(2))*(xi(5)-xi(3))*(xi(5)-xi(4))) * yi(5);
collect(P)

% Polynom mit Vandermonde-Matrix
A = [ones(size(xi)), xi, xi.^2, xi.^3, xi.^4]
b = yi
a = A \ b;
a = flipdim(a,1)    % in poly-Reihenfolge

% Plotwerte
xp1 = 1:0.01:5;
yp1 = polyval(a, xp1);

%% natürliche Splines
A1 = [2 1 0 0 0; 1 4 1 0 0; 0 1 4 1 0; 0 0 1 4 1; 0 0 0 1 2];
b1 = [3 0 -3 9 9]';
M1 = A1 \ b1

spl1 = zeros(4, 4);    % Koeffizienten aller 4 Polynome
for k = 1:4
  s = x - xi(k);
  P = (3 - 2*s)*s^2*yi(k+1) + (1+2*s)*(1-s)^2*yi(k) ...
      + s^2*(s-1)*M1(k+1) + s*(s-1)^2*M1(k);
  P = collect(P);
  spl1(k,:) = sym2poly(P);
end
spl1

% Plotwerte
xp2 = [];
yp2 = [];
for k = 1:4
  xk = k:0.01:(k+1);
  poly = spl1(k,:);
  yk = polyval(poly, xk);
  xp2 = [xp2, xk];
  yp2 = [yp2, yk];
end

%% not-a-knot Splines
A1 = [1 0 -1 0 0; 1 4 1 0 0; 0 1 4 1 0; 0 0 1 4 1; 0 0 1 0 -1];
b1 = [4 0 -3 9 -6]';
M1 = A1 \ b1

spl1 = zeros(4, 4);    % Koeffizienten aller 4 Polynome
for k = 1:4
  s = x - xi(k);
  P = (3 - 2*s)*s^2*yi(k+1) + (1+2*s)*(1-s)^2*yi(k) ...
      + s^2*(s-1)*M1(k+1) + s*(s-1)^2*M1(k);
  P = collect(P);
  spl1(k,:) = sym2poly(P);
end
spl1

% Plotwerte
xp3 = [];
yp3 = [];
for k = 1:4
  xk = k:0.01:(k+1);
  poly = spl1(k,:);
  yk = polyval(poly, xk);
  xp3 = [xp3, xk];
  yp3 = [yp3, yk];
end

% alles plotten
plot(xp1, yp1, 'r-', xp2, yp2, 'c-', xp3, yp3, 'b-', xi, yi, 'k*')
xlim([0.98 5.02]);
legend("Polynom", "natürlicher Spline", "not-a-knot Spline", "Daten", ...
   "Location", "best");
F = getframe(gcf);
imwrite(F.cdata, "bild26.png");

