% Use nature spline to interpolate given data.

% Set n, x_i and f(x_i)
n = 5;
x = [-2, -1, 0, 1, 2, 3];
fx = [1, 4, 11, 16, 13, -4];

% Set the coefficients a_i
a = fx;

% Set h_i. Here we use a Matlab function "diff"
% For more information, type "help diff"
h = diff(x);

% Set the matrix A
A = zeros(n+1,n+1);
A(1,1) = 1;
A(n+1,n+1) = 1;
for i=2:n
    A(i,i-1) = h(i-1);
    A(i,i) = 2*(h(i-1) + h(i));
    A(i,i+1) = h(i);
end

% Set the right-hand side vector rhs
rhs = zeros(n+1,1);
for i=2:n
    rhs(i) = 3*(a(i+1)-a(i))/h(i) - 3*(a(i)-a(i-1))/h(i-1);
end

% Solve for c_i
c = A\rhs;

% solve for b_i and d_i
b = (a(2:(n+1))-a(1:n))./h - h.*(c(2:(n+1))+2*c(1:n))'/3;
d = (c(2:(n+1))-c(1:n))'./h/3;

% Draw the graph
figure(1);
plot(x,fx,'*');
hold on;
for i=1:n
    plotx = x(i):0.1:x(i+1);
    tempx = plotx - x(i);
    ploty = a(i) + b(i)*tempx + c(i)*tempx.^2 + d(i)*tempx.^3;
    plot(plotx, ploty);
end
hold off;