%   The following computes
%     the interpolating polynomial for a given input of data points. Experiment
%     different coordinates (x and y) and try more points. 

 
% set the data points for interpolation
% x-coordinates in x and y-coordinates in y
x=[0; 1; 2;  3;  4];
y=[1; 2; 0; -4; -3];

u=(min(x)-1:0.05:max(x)+1)';  % set points in x-axis

% compute the interpolating polynomial
n = length(x);
v = zeros(size(u));
for k = 1:n
   w = ones(size(u));
   for j = [1:k-1 k+1:n]
      w = (u-x(j))./(x(k)-x(j)).*w;
   end
   v = v + w*y(k);
end

plot(x,y,"+",u,v,"-") 
