% y = polynomial_interpolation(xi, yi, x)
%
% Given a number of sample points listed in a vector xi, and
% associated function values in a same-sized vector yi,
% construct the interpolating polynomial.
% Then, for each point in x, evaluate the polynomial and
% return the result in y.

function y = polynomial_interpolation(xi, yi, x)

% how many points?
n=length(xi);

% form the Vandermonde matrix, A(i,j)=xi(i)^(j-1)
% (using a MATLAB vectorization trick to make this a bit faster - do it row by row)
A=zeros(n,n);
for i=1:n
   A(i,:) = xi(i).^(0:n-1);
end

% get interpolating polynomial coefficients
% (this solves the linear system of inteprolating conditions p(xi)=yi,
% with more MATLAB shortcuts - the backslash operator - so that c contains
% the coefficients of the polynomial in order from constant, to linear, ...,
% up to degree n-1.
c=A\yi(:);

% evaluate the polynomial at requested points
% (this uses more MATLAB vectorization tricks, computing all values one term at a time)
y=zeros(size(x));
for i=1:n
   y=y+c(i)*x.^(i-1);
end
