[proofplan]
We compute the Hamiltonian vector field in the canonical coordinate frame on $T^*Q$. Since the coordinate vector fields form a local frame, we write $X_H$ with unknown coefficient functions in the $\partial_{q_i}$ and $\partial_{p_i}$ directions. Evaluating the defining identity $\omega(X_H,Y)=dH(Y)$ on the coordinate vector fields determines those coefficients. Finally, an integral curve has velocity equal to $X_H$, so the coefficient identities become the coordinate differential equations.
[/proofplan]
[step:Expand the Hamiltonian vector field in the canonical coordinate frame]
On the coordinate domain $V\subset T^*Q$, the coordinate vector fields
$\frac{\partial}{\partial q_1},\dots,\frac{\partial}{\partial q_n},\frac{\partial}{\partial p_1},\dots,\frac{\partial}{\partial p_n}$
form a smooth local frame for $TV$. Therefore there exist smooth coefficient functions $a_i,b_i:V\to\mathbb{R}$ such that
\begin{align*}
X_H=\sum_{i=1}^n a_i\frac{\partial}{\partial q_i}+\sum_{i=1}^n b_i\frac{\partial}{\partial p_i}.
\end{align*}
Define the canonical symplectic form $\omega$ on $V$ by
\begin{align*}
\omega=\sum_{i=1}^n dq_i\wedge dp_i.
\end{align*}
With the convention used here, the Hamiltonian vector field $X_H$ is the unique vector field satisfying
\begin{align*}
\omega(X_H,Y)=dH(Y)
\end{align*}
for every smooth vector field $Y$ on $V$. We will determine $a_i$ and $b_i$ from this defining equation for $X_H$.
[/step]
[step:Evaluate the defining equation on $\frac{\partial}{\partial p_j}$ to identify the $q_j$ coefficient]
Fix $j\in\{1,\dots,n\}$. Since the coordinate one-forms are dual to the coordinate vector fields, for all $i,k\in\{1,\dots,n\}$ we have
\begin{align*}
dq_i\left(\frac{\partial}{\partial q_k}\right)=\delta_{ik}.
\end{align*}
Also
\begin{align*}
dq_i\left(\frac{\partial}{\partial p_k}\right)=0.
\end{align*}
Similarly,
\begin{align*}
dp_i\left(\frac{\partial}{\partial q_k}\right)=0.
\end{align*}
Finally,
\begin{align*}
dp_i\left(\frac{\partial}{\partial p_k}\right)=\delta_{ik}.
\end{align*}
Using these identities, we compute
\begin{align*}
(dq_i\wedge dp_i)\left(X_H,\frac{\partial}{\partial p_j}\right)=dq_i(X_H)\,dp_i\left(\frac{\partial}{\partial p_j}\right)-dq_i\left(\frac{\partial}{\partial p_j}\right)\,dp_i(X_H).
\end{align*}
The first factor gives $dq_i(X_H)=a_i$, the second gives $dp_i(\frac{\partial}{\partial p_j})=\delta_{ij}$, and the term $dq_i(\frac{\partial}{\partial p_j})$ is zero. Hence
\begin{align*}
(dq_i\wedge dp_i)\left(X_H,\frac{\partial}{\partial p_j}\right)=a_i\delta_{ij}.
\end{align*}
Summing over $i$ gives
\begin{align*}
\omega\left(X_H,\frac{\partial}{\partial p_j}\right)=a_j.
\end{align*}
By the defining identity $\omega(X_H,Y)=dH(Y)$, applied to $Y=\frac{\partial}{\partial p_j}$, we also have
\begin{align*}
\omega\left(X_H,\frac{\partial}{\partial p_j}\right)=dH\left(\frac{\partial}{\partial p_j}\right)=\frac{\partial H}{\partial p_j}.
\end{align*}
Therefore
\begin{align*}
a_j=\frac{\partial H}{\partial p_j}.
\end{align*}
[guided]
Fix $j\in\{1,\dots,n\}$. The purpose of testing against $\frac{\partial}{\partial p_j}$ is to isolate the coefficient of $X_H$ in the $\frac{\partial}{\partial q_j}$ direction. The coordinate one-forms act on the coordinate vector fields by the duality identities, valid for all $i,k\in\{1,\dots,n\}$:
\begin{align*}
dq_i\left(\frac{\partial}{\partial q_k}\right)=\delta_{ik}.
\end{align*}
Also
\begin{align*}
dq_i\left(\frac{\partial}{\partial p_k}\right)=0.
\end{align*}
Similarly,
\begin{align*}
dp_i\left(\frac{\partial}{\partial q_k}\right)=0.
\end{align*}
Finally,
\begin{align*}
dp_i\left(\frac{\partial}{\partial p_k}\right)=\delta_{ik}.
\end{align*}
Using the definition of the wedge product of one-forms, for each $i$ we have
\begin{align*}
(dq_i\wedge dp_i)\left(X_H,\frac{\partial}{\partial p_j}\right)=dq_i(X_H)\,dp_i\left(\frac{\partial}{\partial p_j}\right)-dq_i\left(\frac{\partial}{\partial p_j}\right)\,dp_i(X_H).
\end{align*}
From the expansion
\begin{align*}
X_H=\sum_{k=1}^n a_k\frac{\partial}{\partial q_k}+\sum_{k=1}^n b_k\frac{\partial}{\partial p_k},
\end{align*}
we get $dq_i(X_H)=a_i$. Also $dp_i(\frac{\partial}{\partial p_j})=\delta_{ij}$ and $dq_i(\frac{\partial}{\partial p_j})=0$. Thus
\begin{align*}
(dq_i\wedge dp_i)\left(X_H,\frac{\partial}{\partial p_j}\right)=a_i\delta_{ij}.
\end{align*}
Now sum this identity over $i$ because the canonical symplectic form was defined by $\omega=\sum_{i=1}^n dq_i\wedge dp_i$. Only the term $i=j$ survives, so
\begin{align*}
\omega\left(X_H,\frac{\partial}{\partial p_j}\right)=a_j.
\end{align*}
With the convention used in this theorem, the Hamiltonian vector field is defined by the identity $\omega(X_H,Y)=dH(Y)$ for every smooth vector field $Y$ on $V$. Applying this identity with $Y=\frac{\partial}{\partial p_j}$ gives
\begin{align*}
a_j=\omega\left(X_H,\frac{\partial}{\partial p_j}\right)=dH\left(\frac{\partial}{\partial p_j}\right)=\frac{\partial H}{\partial p_j}.
\end{align*}
So the $\frac{\partial}{\partial q_j}$ coefficient of $X_H$ is exactly the partial derivative of $H$ with respect to $p_j$.
[/guided]
[/step]
[step:Evaluate the defining equation on $\frac{\partial}{\partial q_j}$ to identify the $p_j$ coefficient]
Fix $j\in\{1,\dots,n\}$. For each $i$,
\begin{align*}
(dq_i\wedge dp_i)\left(X_H,\frac{\partial}{\partial q_j}\right)=dq_i(X_H)\,dp_i\left(\frac{\partial}{\partial q_j}\right)-dq_i\left(\frac{\partial}{\partial q_j}\right)\,dp_i(X_H).
\end{align*}
Here $dq_i(X_H)=a_i$, $dp_i(\frac{\partial}{\partial q_j})=0$, $dq_i(\frac{\partial}{\partial q_j})=\delta_{ij}$, and $dp_i(X_H)=b_i$. Hence
\begin{align*}
(dq_i\wedge dp_i)\left(X_H,\frac{\partial}{\partial q_j}\right)=-\delta_{ij}b_i.
\end{align*}
Summing over $i$ gives
\begin{align*}
\omega\left(X_H,\frac{\partial}{\partial q_j}\right)=-b_j.
\end{align*}
By the defining identity for $X_H$,
\begin{align*}
-b_j=dH\left(\frac{\partial}{\partial q_j}\right)=\frac{\partial H}{\partial q_j}.
\end{align*}
Therefore
\begin{align*}
b_j=-\frac{\partial H}{\partial q_j}.
\end{align*}
[/step]
[step:Translate the vector field identity into differential equations for integral curves]
Combining the coefficient identities for all $i\in\{1,\dots,n\}$, we obtain
\begin{align*}
X_H=\sum_{i=1}^n \frac{\partial H}{\partial p_i}\frac{\partial}{\partial q_i}-\sum_{i=1}^n \frac{\partial H}{\partial q_i}\frac{\partial}{\partial p_i}.
\end{align*}
Let $\gamma:I\to V$ be an integral curve of $X_H$. By definition of integral curve,
\begin{align*}
\dot\gamma(t)=X_H(\gamma(t))
\end{align*}
for every $t\in I$. In canonical coordinates,
\begin{align*}
\dot\gamma(t)=\sum_{i=1}^n \dot q_i(t)\frac{\partial}{\partial q_i}\bigg|_{\gamma(t)}+\sum_{i=1}^n \dot p_i(t)\frac{\partial}{\partial p_i}\bigg|_{\gamma(t)}.
\end{align*}
Evaluating the formula for $X_H$ at $\gamma(t)$ and comparing coefficients in the coordinate frame gives, for every $i\in\{1,\dots,n\}$,
\begin{align*}
\dot q_i(t)=\frac{\partial H}{\partial p_i}(\gamma(t)).
\end{align*}
and
\begin{align*}
\dot p_i(t)=-\frac{\partial H}{\partial q_i}(\gamma(t)).
\end{align*}
These are Hamilton's equations in the canonical coordinates induced by $(q_1,\dots,q_n)$.
[/step]
