Interpolating polynomial, Vandermonde method

intpoly1.m

x = [1 2 3 4]'; % transpose to make a column 
y = [3 1 4 1]'; 
n = length(x);  % number of elements 
M = zeros(n,n); % initialize a matrix n by n
for k=1:n
    M(:,k) = x.^(k-1);    % fill the matrix with powers of x-values
end
coeff = M\y;    % find the coefficients 
t = linspace(1,4);
poly = coeff(n)*ones(size(t));  % prepare for evaluation 
for k=1:n-1
    poly=poly.*t + coeff(n-k);  % nested multiplication
end
plot(t, poly)
hold on
plot(x, y, 'r*')          % add original points to show the curve passes through them