function y=lagrange(xdat,ydat,x)
% Computes the Lagrange polynomial
% using the datapoints in (xdat, ydat)

N=length(x);
M=length(xdat);
D=zeros(M,M);
y=zeros(1,N);
for k=1:N
    xtst=x(k);
    [xb,ind]=sort(abs(xtst-xdat));
    xsort=xdat(ind);
    D(:,1)=ydat(ind)';
    for j=1:M
        for i=1:M-j
            D(i,j+1)=(D(i+1,j)-D(i,j))/(xsort(i+j)-xsort(i));
        end
    end
    xdiff=xtst-xsort;
    prod=1;
    for i=1:M
        y(k)=y(k)+prod*D(1,i);
        prod=prod*xdiff(i);
    end
end