n = 10;
x = [-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1];
fx = [0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385];
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.01: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;
axis([-1, 1, -0.5,2]);