function [d,z] = div_dif_hermite(x,y,ypr)
%
% Given vectors of interpolation points x and interpolation values 
% y and ypr, this produces the vector of divided differences d and 
% the vector of points z used in constructing the Newton form of 
% the Hermite interpolating polynomial H(x) satisfying 
% 
%   H(x(i)) = y(i),  H'(x(i)) = ypr(i)  for i = 1,...,n, 
%
% where n = length(x) = length(y) = length(ypr). The form of H(x) is 
%
% H(x) = d(1) + d(2)*(x-z(1)) + ... + d(2*n)*(x-z(1))*...(x-z(2*n-1)) 
% 
% To evaluate H(x) for a given x, first use this routine to compute 
% d and z and then compute H(x) = Hval with 
%
%                 Hval = p_eval(x,z,d); 
%
% where p_eval is the routine in p_eval.m. 
%
if size(x) ~= size(y)|size(x) ~= size(ypr)
  error('The input vectors x, y, and ypr must have the same size.'); 
end
n = length(x);
nn = 2*n; 

% Row or column vectors? 

if n == size(x,2)
  row = 1;
else
  row = 0;
end

% Create z. 

if row
  z = reshape([x;x],1,nn);
else
  z = reshape([x';x'],nn,1);
end 

% Create d with d(1) = y(1) and the first-order differences in d(2:n). 

y(2:n) = (y(2:n) - y(1:n-1))./(x(2:n) - x(1:n-1)); 

if row 
  d = reshape([y;ypr],1,nn);
else
  d = reshape([y';ypr'],nn,1);
end 

% Create the remaining differences in d(3:n). 
  
for j = 3:nn 
  d(j:nn) = (d(j:nn) - d(j-1:nn-1))./(z(j:nn) - z(1:nn-j+1)); 
end
