function y=hermfun(X,Y,YP,x)
H=hermite(X,Y,YP);
n=length(X);
y=0;
for i=1:n
    p1=H(2*i-1,2*i-1);
    p2=H(2*i,2*i).*(x-X(i));
    for j=1:i-1
        p1=p1.*(x-X(j)).^2;
        p2=p2.*(x-X(j)).^2;
    end
    y=y+p1+p2;
end
       
% hermfun([1.3 1.6 1.9],[0.6200860 0.4554022 0.2818186],[-0.5220232 -0.5698959 -0.5811571],1.5)