for i = 1:n,
 L(i) = 1
end

i = 0
for j = i+1:n
    L(i) = L(i)*(x-x(j))/(x(i)-x(j))
end

for i = 1:n-1
  for j = 0:i-1
   L(i) = L(i)*(x-x(j))/(x(i)-x(j))
  end
  for j = i+1:n
   L(i) = L(i)*(x-x(j))/(x(i)-x(j))
  end
end

i = n
for j = 0:i-1
    L(i) = L(i)*(x-x(j))/(x(i)-x(j))
end

% the above generates L_i,n(x)
% next compute P(x)

P = 0
for i = 0:n
 P = P + f(i)*L(i)
end