function F=hermite_table(xi,fi,dfi)

% Hermite interpolation (variation on Newton's divided difference): See section 3.3.
% This version returns the table of computed values, as in Algorithm 3.3.
% First input is the vector of nodes (sample x-coordinates).
% Second input is either the corresponding vector of y-coordinates,
%  or the function f.
% Third input is either the corresponding vector of values of y',
%  or the function f'.

n=length(xi)-1;

% Create the vector z = [x_0 x_0 .. x_n x_n].
zi(1:2:2*n+1) = xi(:);
zi(2:2:2*n+2) = xi(:);

% Fill the first column of the table F using y-data or the function f.
if isa(fi, 'function_handle')
    F(1:2*n+2, 1) = fi(zi(:));
else
    F(1:2:2*n+1, 1) = fi(:);
    F(2:2:2*n+2, 1) = fi(:);
end
% Fill half of the 2nd column of the table F using y'-data or the function f'.
if isa(fi, 'function_handle')
    F(2:2:2*n+2, 2) = dfi(xi(:));
else
    F(2:2:2*n+2, 2) = dfi(:);
end

for i=1:2*n+1
    for j=1:i
        if j > 1  |  mod(i, 2) == 0 %If j>1 or i is even
            F(i+1, j+1) = (F(i+1,j) - F(i,j)) ...
                         /(zi(i  +1)-zi(i-j+1));
        end
    end
end