function f = HermInt(x, y, dy)
%**************************************************************************
% f = HermInt(x, y, dy)
%
% The Hermite polynomial interpolant (first derivatives are given).
% x is a row n-vector with distinct components, y is a matrix of row
% n-vectors evaluated at the x components, and dy is a matrix of first
% derivative values computed at the corresponding y-values.  f is a row
% n-vector with the property that if
%
%      p(x) = f(1) + f(2)(x-x(1)) + f(3)(x-x(1))^2 + f(4)(x-x(1))^2*
%             *(x-x(2) + ... + f(2*n)(x-x(1))^2*...*(x-x(n-2))^2*(x-x(n-1))
% then
%      p(x(2*i)) = y(2*i), i=1:n.
%**************************************************************************

[m, n] = size(y); [m2, n2] = size(dy); [one n3] = size(x);
if (m - m2) | (n - n2) | (n - n3) | (one - 1)   % Are sizes the same?
    f = NaN;                                    % If not, assign NaN to "c"

else
    x = reshape([x; x], 1, 2*n);        % Cause x to repeat number pairs
    n = 2*n;                            % Double the number of points
    
    f = zeros(m, n);    % Start with empty matrix
    f(:,1) = y(:,1);    % Make first column the same as that of y

    % This first FOR loop creates a matrix c with column vectors that
    % alternate between y(i)' and y[i,i+1] (and c(:,1) = y(:,1) above)
    for p = 2:2:(n-2)
           f(:,p)   = dy(:,p/2);
           f(:,p+1) = (y(:,p/2+1) - y(:,p/2)) ./ (x(p+1) - x(p));
    end
    f(:,n) = dy(:,n/2);

    % This next set of FOR loops computes further divided differences
    % normally (InterpN.m does these when x and y are column vectors).
    for k = 2:(n-1)
        for p = n:-1:(k+1)
            f(:,p) = (f(:,p) - f(:,p-1)) ./ (x(p) - x(p-k));
        end
    end
end