function p = polyinterp(xdata, fdata, x)
% POLYINTERP interpolates the given data using the
%            Lagrange Interpolation formula
%
% xdata - the x coordinates of the nodes 
% fdata - the function value at the nodes
% x     - 'fine' mesh
% p     - interpolated data on the 'fine' mesh
%
	M = length(xdata);     
	s = size(x);

    p = zeros(s);
	for k = 1:M
	    w = ones(s);
	    for j = [1:k-1 k+1:M]
			w = (x - xdata(j)) ./ (xdata(k) - xdata(j)) .* w;
	    end
	    p = p + w*fdata(k);
	end
end
