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

function H = HR0_U(U, F)
A = U + F;
[n, n] = size(A);
In = eye(n); 
In2 = eye(n^2);

F = A-U;
N = inv(In - U);
R = F*N;
Rt = R';
[R0, wr, vr] = domeig(R);

DRU = kron(N', R); % Eq (50)
HR0R = Hlambda_A(R); % H[R0,vecR], using Eq (24)

HRU = kron(kron(In, Kmn(n,n)), In) * (kron(In2,Rt(:))*kron(N',N) + kron(N(:),In2)*Kmn(n,n)*kron(N', R)); % Eq (56)

B = kron(kron(wr', vr'), In2)*HRU + DRU'*HR0R*DRU; % Eq (57)
H = 1/2*(B + B');

end