function S = cubicsplineb(t,y)
%CUBICSPLINEB  Natural cubic spline interpolation.
%
%   S = cubicsplineb(t,y) computes a natural cubic spline approximation
%   to the data t, y. It then plots the data and the cubic spline.
%   (ti, yi) contains n+1 data points. In the n x 4 output matrix S,
%   row 1 is the cubic polynomial fit to t1 <= t <= t2, ...,
%   Row n is the cubic polynomial fit to tn <= t <= tnp1.
%   The function makes a plot of the data and the cubic spline
%   approximation. 

% Build the diagonals of the (n+1) x (n+1) tridiagonal matrix.
n = length(t)-1;
a = ones(n,1);
b = 4*ones(n+1,1);
c = a;
b(1) = 2;
b(n+1) = 2;

%   Build the right-hand side of the system.
rhs = zeros(n+1,1);
rhs(1) = 3*(y(2) - y(1));
rhs(n+1) = 3*(y(n+1) - y(n));
for i = 2:n
   rhs(i) = 3*(y(i+1) - y(i-1));
end

% Solve the system using the linear Thomas algorithm
D = thomas(a,b,c,rhs);

% onstruct the cubic polynomials in the rows of S using
% MATLAB polynomial format.
S = zeros(n,4);
for i = 1:n
      S(i,1) = 2*(y(i) - y(i+1)) + D(i) + D(i+1);
      S(i,2) = 3*(y(i+1) - y(i)) - 2*D(i) -D(i+1);
      S(i,3) = D(i);
      S(i,4) = y(i);
end

% Plot the data and the cubic spline approximation.
hold on;
%   Each cubic is parameterized over 0 <= u <= 1.
u = 0.0:.01:1.0;
for i = 1:n
   p = S(i,:);
   % corresponding t values.
   x = t(i):(t(i+1)-t(i))/100:t(i+1);
   plot(x,polyval(p,u),'k');
end

plot(t,y,'*r');
hold off;
