
% Cubic spline interpolation --  with n=3 as we did in class

% set up the data points (x(j),y(j)) where x and y are defined below 
x=[0 1 3];
y=[0 3 1]

%%% define the matrix S as in Section I.2.6 (p. 42 in Ch 1 of the notes)
% initially build the 8 x 8 zero matrix
S=zeros(8,8);
% then define the non-zero entries for the first two rows
S(1,1:4)=x(1).^(3:-1:0);
S(2,1:4)=x(2).^(3:-1:0);

% Alternatively
for j=1:2; 
    S(j,1:4)=x(j).^(3:-1:0);
end

% next, rows 3 and 4
for j=3:4;
S(j,5:end)=x(j-1).^(3:-1:0);
end

% finally, the remaining rows
S(5,:)=[3*x(2)^2 2*x(2) 1 0 -3*x(2)^2 -2*x(2) -1 0]
S(6,:)=[6*x(2) 2 0 0 -6*x(2) -2 0 0]
S(7,:)=[6*x(1) 2 0 0 0 0 0 0]
S(8,:)=[0 0 0 0 6*x(3) 2 0 0]


% define the vector b 
b=[y(1); y(2); y(2); y(3); 0; 0; 0; 0];

% solve for a where a(1)=A_1, a(2)=A_2, ...
a=S\b;

% plot the resulting cubic spline in a new figure;
figure;
plot(x,y, '*','MarkerSize',20,'LineWidth',2);
hold on
X1=linspace(x(1),x(2),100);
plot(X1, polyval(a(1:4),X1),'r','LineWidth',2);
X2=linspace(x(2),x(3),100);
plot(X2, polyval(a(5:8),X2),'b','LineWidth',2);

% check the condition number of S
cond(S)

%% Now repeat after shifting the x coordinates of the data by 100


% set up the data points (x(1),y(1)) where x and y are defined below 
x=[0 1 3]+100;
y=[0 3 1]

%%% define the matrix S as in Section I.2.6 (p. 42 in Ch 1 of the notes)
% initially build the 8 x 8 zero matrix
S=zeros(8,8);
% then define the non-zero entries for the first two rows
S(1,1:4)=x(1).^(3:-1:0);
S(2,1:4)=x(2).^(3:-1:0);

% Alternatively
for j=1:2; 
    S(j,1:4)=x(j).^(3:-1:0);
end

% next, rows 3 and 4
for j=3:4;
S(j,5:end)=x(j-1).^(3:-1:0);
end

% finally, the remaining rows
S(5,:)=[3*x(2)^2 2*x(2) 1 0 -3*x(2)^2 -2*x(2) -1 0]
S(6,:)=[6*x(2) 2 0 0 -6*x(2) -2 0 0]
S(7,:)=[6*x(1) 2 0 0 0 0 0 0]
S(8,:)=[0 0 0 0 6*x(3) 2 0 0]


% define the vector b 
b=[y(1); y(2); y(2); y(3); 0; 0; 0; 0];

% solve for a where a(1)=A_1, a(2)=A_2, ...
a=S\b;

% plot the resulting cubic spline in a new figure;
figure;
plot(x,y, 'o','MarkerSize',20,'LineWidth',2);
hold on
X1=linspace(x(1),x(2),100);
plot(X1, polyval(a(1:4),X1),'LineWidth',2);
X2=linspace(x(2),x(3),100);
plot(X2, polyval(a(5:8),X2),'LineWidth',2);

% check the condition number of S
cond(S)


