% Cubic Spline Interpolation - Clamped
% INPUT: xi is the grid (points on x-axis) and a are points on y-axis. 
% fp is a two vector [f'(x0),f'(xN)] specifying the derivative at each end.
 
function [P] = cubic_splineclp(xi,a,fp)
 
if length(xi) ~= length(a)
 erro('vectors xi and a must be of same length');
end
 
% Plotting points we want to interpolate between:
grid on; hold on;
title('Cubic Spline Interpolation');
plot(xi,a,'or');
 
n = length(xi);
 
% Vector h with subintervals:
h = zeros(n-1,1);
for j = 1:n-1
 h(j) = xi(j+1) - xi(j);
end
 
% Coefficient matrix A:
A = zeros(n);
 
% Clamped boundary conditions:
A(1,1) = 2*h(1);
A(1,2) = h(1);
A(n,n) = 2*h(n-1);
A(n,n-1) = h(n-1);
 
for i = 2:n-1
 A(i,i-1) = h(i-1);
 A(i,i) = 2*(h(i-1)+h(i));
 A(i,i+1) = h(i);
end
 
% Vector b:
b = zeros(n,1);
b(1) = (3/h(1))*(a(2)-a(1)) - 3*fp(1);
b(n) = 3*fp(2) - (3/h(n-1))*(a(n)-a(n-1));
 
for i = 2:n-1
 b(i) = (3/h(i))*(a(i+1)-a(i)) - (3/h(i-1))*(a(i)-a(i-1));
end
 
% Coefficient vector cj:
cj = A\b;
 
% Coefficient vector bj:
bj = zeros(n-1,1);
for i = 1:n-1
 bj(i) = (1/h(i))*(a(i+1)-a(i)) - (1/3*h(i))*(2*cj(i)+cj(i+1));
end
 
% Coefficient vector dj:
dj = zeros(n-1,1);
for i = 1:n-1
 dj(i) = (1/(3*h(i))) * (cj(i+1)-cj(i));
end
 
% Making a matrix P with all polynomials (and intervals)
P = zeros(n-1,4);
X = zeros(n-1,2);
for i = 1:n-1
 P(i,1) = dj(i);
 P(i,2) = cj(i);
 P(i,3) = bj(i);
 P(i,4) = a(i);
 X(i,1) = xi(i);
 X(i,2) = xi(i+1);
end
 P = [P,X];
 
% Plotting results:
resolution = 100;
for i = 1:n-1
 f = @(x) a(i) + bj(i).*(x-xi(i)) + cj(i).*(x-xi(i)).^2 + dj(i).*(x-xi(i)).^3;
 xf = linspace(xi(i),xi(i+1),resolution);
 plot(xf,f(xf),'m-');
end
print -depsc spline_clp.eps
end