% program developed by M. R. Grasselli & D. E. Pelinovsky
% as part of the courseware "Numerical Interpolation &
% Approximation Theory" (2006)


close all;
clear all;

x = [0 0.2 0.4 0.6 0.8 1.0];
y = [2.0079    2.3454    3.0197    4.1274    5.7655    8.0315];
xInt = 0 : 0.01 : 1;
% degree of the interpolating polynomial
n = length(x) - 1; 
ni = length(xInt);
% First step: prepare a table of Lagrange polynomials
L = ones(ni,n+1); 
for j = 1 : (n+1) 
    for i = 1 : (n+1) 
        if (i ~= j) 
            L(:,j) = L(:,j).*(xInt' - x(i))/(x(j)-x(i));
        end
    end
end
% Second step: compute the interpolating polynomial at the set of x-values
yInt = y*L'; 
plot(xInt,yInt,'g',x,y,'b*');
