n = 5;
x = [-2, -1, 0, 1, 2, 3];
fx = [1, 4, 11, 16, 13, -4];
a = fx;
h = diff(x);
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
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
c = A\rhs;
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;
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;