function c = hermdiff(x,y,slope)
%  This function returns the coefficients for the
%  Newton's divided difference form of the
%  Hermite polynomial
%  Input:  vectors x and y of data points
%             vector slope of slope values

m = max(size(x));
deg = 2*m-1;

%  First, stick in the y-values
for k=1:m
  c(2*k-1)= y(k);
  c(2*k) = y(k);
  xx(2*k-1)= x(k);
  xx(2*k) = x(k);
end

%  Next, deal with the first divided differences only computing every other term.

for j = deg:-2:3
    c(j) = (c(j)-c(j-1))/(xx(j)-xx(j-1));
end

c

%  Then, fix those first divided differences which should be slopes
for j=2:2:deg+1
  c(j) = slope(j/2);
end

c

%  Finally, do the rest of the divided difference table as before.

for k=2:deg
  for j=deg+1:-1:k+1
    c(j) = (c(j)-c(j-1))/(xx(j)-xx(j-k));
  end
end