function Pz=lagrange(x,y,z)
% Evaluates the nth Lagrange interpolating polynomial P 
%     to the function f using (n+1) nodes 
% Input x: row vector containing (n+1) nodes
% Input y: row vector containing function values at the nodes
% Input z: point you are interested in approximating the function at
% Output: value of P(z),  

Pz=0;
n=length(x);
L=ones(1,n);

L(1) = prod(z-[x(2:n)]) / ( prod(x(1) -[x(2:n) ]));

L(n) = prod(z-[x(1:n-1)]) / ( prod(x(n) -[x(1:n-1) ]));

for i = [2:n-1]
L(i) = prod(z-[x(1:i-1)]).*prod(z-[x(i+1:n)])/(prod(x(i)-[x(1:i-1)]).*prod(x(i)-[x(i+1:n)]));
end

Pz=sum(L.*y);

end