function y=hmt(n,x,f,f1,xi)
y=0.0;
for i=1:n
    h=1.0;
    a=0.0;
    for j=1:n
        if j~=i
            h=h*((xi-x(j))/(x(i)-x(j)))^2;
            a=1/(x(i)-x(j))+a;
        end
    end
    y=y+h*((x(i)-xi)*(2*a*f(i)-f1(i))+f(i));
end