function d=hermitedif(x,y,yp)

%tic

n=length(x);
nn=2*n;

for i=1:n
  xx(2*i-1)=x(i);
  xx(2*i)=x(i);
end

d(1)=y(1);

for j=1:n
  d(2*j)=yp(j);
end

for j=n:-1:2
  d(2*j-1)=(y(j)-y(j-1))/(x(j)-x(j-1));;
end

for i=2:nn-1
   for j=nn:-1:i+1
      d(j)=(d(j)-d(j-1))/(xx(j)-xx(j-i));  
   end
end
    
%toc