  We consider a scalar function depending on a numerical solution of an initial
value problem, and its second-derivative (Hessian) matrix for the initial
value. The need to extract the information of the Hessian or to solve a linear
system having the Hessian as a coefficient matrix arises in many research
fields such as optimization, Bayesian estimation, and uncertainty
quantification. From the perspective of memory efficiency, these tasks often
employ a Krylov subspace method that does not need to hold the Hessian matrix
explicitly and only requires computing the multiplication of the Hessian and a
given vector.
  One of the ways to obtain an approximation of such Hessian-vector
multiplication is to integrate the so-called second-order adjoint system
numerically. However, the error in the approximation could be significant even
if the numerical integration to the second-order adjoint system is sufficiently
accurate. This paper presents a novel algorithm that computes the intended
Hessian-vector multiplication exactly and efficiently. For this aim, we give a
new concise derivation of the second-order adjoint system and show that the
intended multiplication can be computed exactly by applying a particular
numerical method to the second-order adjoint system. In the discussion,
symplectic partitioned Runge--Kutta methods play an essential role.
