I am reading a book about geometric numerical integration. In the text, they present Euler symplectic integration
\begin{equation} p_{n+1}=p_n-hH_q(p_{n+1},q_n)\\ q_{n+1} = q_n+hH_p(p_{n+1},q_n) \end{equation}
for integrating the Hamiltonian system \begin{align} \dot{p}=-H_q(p,q)\\ \dot{q}=H_p(p,q) \end{align}
In the book it is written that differentiating the first system (symplectic integration scheme) with respect to $(p_n,q_n)$ yields
$$\begin{bmatrix} I+hH_{qp}^T &0\\ -hH_{pp}& I \end{bmatrix} (\frac{\partial(p_{n+1},q_{n+1})}{\partial (p_n,q_n)}) = \begin{bmatrix} I &-H_{qq}\\ 0& I+hH_{qp} \end{bmatrix}$$
I think I understand where the matrix on the right side comes from since it is the jacobian matrix of $$f(p_n,q_n) =\begin{bmatrix}p_n-hH_q(p_{n+1},q_n)\\ q_n+hH_p(p_{n+1},q_n)\end{bmatrix}$$ with respect to $(p_n,q_n)$ but I don't really understand how we get the expression on the left, i.e. where does
$$\begin{bmatrix} I+hH_{qp}^T &0\\ -hH_{pp}& I \end{bmatrix} (\frac{\partial(p_{n+1},q_{n+1})}{\partial (p_n,q_n)}) $$ come from?
You have to fully write out the total differentials of the equations, then you get \begin{align} dp_{n+1}&=dp_n-hH_{qp}dp_{n+1}-hH_{qq}dq_n\\ dq_{n+1}&=dq_n+hH_{pp}dp_{n+1}+hH_{pq}dq_n \end{align} Now sort the components with the same index to each side and write it in matrix form $$ \pmatrix{I+hH_{qp}&0\\-hH_{pp}&I}\pmatrix{dp_{n+1}\\dq_{n+1}} = \pmatrix{I&-hH_{qq}\\0&I+H_{pq}}\pmatrix{dp_n\\dq_n} $$ Then insert the vector field $(\frac\partial{\partial p_n},\frac\partial{\partial q_n})$ to get the claimed identity.