function poly=hermite(X,Y,M)
%Y is f(X),M is f'(X). if not given then the value is nan
%length(X)=length(Y)=length(M)
n=sum(~isnan(Y));
m=sum(~isnan(M));
p=n+m;
A=zeros(p,p);
b=zeros(p,1);
row=1;
for h=1:length(X)
    if ~isnan(Y(h))
      A(row,:)=X(h).^(0:p-1);
      b(row)=Y(h);
      row=row+1;
    end;  
    if ~isnan(M(h))
        % a0*1+a1*x^1+a2*x^2+a3*x^3=[1 x,x^2,X^3]*[a0;a1;a2;a3]
        A(row,:)=[0,(X(h).^(0:p-2)).*(1:p-1)];
        b(row)=M(h);
        row=row+1;
    end;    
end;
poly=A\b;
return;