% y = lagrange(wk,xk,x)
% arguments:
%   wk : row vector of size N+1 containing the interpolation coefficients
%   xk : row vector of size N+1 containing the interpolation points, xk(1)<...<xk(N+1)
%   x  : points at which the Lagrange polynomial is to be evaluated
% output:
%   y  : value of the Lagrange polynomial at x

function y = lagrange(wk,xk,x)

    [dum,N] = size(wk);
    y = 0.;
    for k = 1:N
        prod = wk(k);
        for j = 1:N
            if j ~= k
                prod = prod * (x-xk(j)) / (xk(k)-xk(j));
            end
        end
        y = y + prod;
    end
  
end