function result = poly_interp(x, y)
% x and y are column vectors with the x and y values of the data points

% there are n+1 data points
n = length(x) - 1;

% construct the Vandermonde matrix
V = zeros(n+1,n+1);
for i=1:n+1
    for j=1:n+1
        V(i,j) = x(i).^(j-1);
    end %for
end %for

V

% solve the system of equations
alpha = V\y

%reverse the alpha vector to match Matlab standards 
%for polynomial coefficient vectors
result = fliplr(alpha');

% plot the solution
xx = x(1):0.01:x(n+1);
yp = polyval(result,xx);

figure(100)
hold off
plot(x,y,'bo')
hold on
plot(xx,yp,'r')
hold off

end %function