function ShowCubicInterp()
% Polynomial interpolation as a linear equation problem.
%
% Given (x(i),y(i)), i=1:4, find a(1:4) so that if
% p(z) = a(1) + a(2)z + a(3)z^2 + a(4)z^3
% then p(x(i)) = y(i), i=1:4.
clc
close all
% Sample problem...
x = [-1 1 2 4]';
y = [ 4 -2 1 -2]';

% Must have
% y(1) = p(x(1)) = a(1) + a(2)x(1) + a(3)x(1)^2 + a(4)x(1)^3
% y(2) = p(x(2)) = a(1) + a(2)x(2) + a(3)x(2)^2 + a(4)x(2)^3
% y(3) = p(x(3)) = a(1) + a(2)x(3) + a(3)x(3)^2 + a(4)x(3)^3
% y(4) = p(x(4)) = a(1) + a(2)x(4) + a(3)x(4)^2 + a(4)x(4)^3
% i.e., y = V*a where
V = [ 1 x(1) x(1)^2 x(1)^3;...
1 x(2) x(2)^2 x(2)^3;...
1 x(3) x(3)^2 x(3)^3;...
1 x(4) x(4)^2 x(4)^3];
a = V\y
xVals = linspace(-1.1,4.1)';
pVals = ( ( a(4)*xVals + a(3) ).*xVals + a(2) ).*xVals + a(1);
plot(xVals,pVals,x,y,'or',[-2 5],[0 0],'-m')
shg

function V = MakeV(x)
% x is a column n-vector
% V is the nxn matrix where V(i,j) = x(i)^(j-1)
n = length(x);
V = zeros(n,n);
V(:,1) = ones(n,1);
for j=2:n
V(:,j) = x.*V(:,j-1);
end