% Function to calculate the Hessian of R0 to the entries of fertility
% matrix F.
%
% Input:
%    U = the transition matrix (n by n matrix)
%    F = the fertility matrix (n by n matrix)
%
% Function outputs result of Equation (63) (n^2 by n^2 matrix).
%

function H = HR0_F(U, F)
A = U + F;
[n, n] = size(A);
In = eye(n); 

F = A-U;
N = inv(In - U);
R = F*N;
HR0R = Hlambda_A(R); % H[R0,vecR], using Eq (24)

B = kron(N,In) * HR0R * kron(N',In); % Eq (62)
H = 1/2*(B + B');

end