function [p] = lip(xp,yp,x)
% Lagrange interploating polynomial
% xp = [x0,x1,...,xn]  vector 0f nodes
% yp = [y0,y1,...,yn]  vector of function values
% x = point at which to evaluate the polynomial

n = length(xp) - 1;
p = 0;

for k = 0 : n
   term = yp(k+1);
   for j = 0 : n
      if j ~= k
         term = term * (x - xp(j+1)) / (xp(k+1) - xp(j+1));
      end
   end
   p = p + term;
end