[step:Compute the bracket of product left-invariant vector fields componentwise]
Let $X_1,X_2\in\mathfrak g$ and $Y_1,Y_2\in\mathfrak h$. Let $\widetilde X_1,\widetilde X_2\in\mathfrak X(G)$ and $\widetilde Y_1,\widetilde Y_2\in\mathfrak X(H)$ be their corresponding left-invariant vector fields.
Define vector fields $A,B\in\mathfrak X(K)$ by
\begin{align*}
A_{(g,h)}=(\widetilde X_{1,g},\widetilde Y_{1,h})
\end{align*}
and
\begin{align*}
B_{(g,h)}=(\widetilde X_{2,g},\widetilde Y_{2,h})
\end{align*}
for every $(g,h)\in K$. By the previous step, $A=\widetilde{(X_1,Y_1)}$ and $B=\widetilde{(X_2,Y_2)}$.
Choose a product coordinate chart around $(g,h)$ with coordinates $(x_1,\dots,x_m,y_1,\dots,y_n)$, where $(x_1,\dots,x_m)$ are coordinates from a chart on $G$ and $(y_1,\dots,y_n)$ are coordinates from a chart on $H$. In these coordinates, there are smooth coefficient functions $a_i,b_i$ on the $G$-chart and $\alpha_j,\beta_j$ on the $H$-chart such that
\begin{align*}
A = \sum_{i=1}^{m} a_i(x)\,\partial_{x_i} + \sum_{j=1}^{n} \alpha_j(y)\,\partial_{y_j}
\end{align*}
and
\begin{align*}
B = \sum_{i=1}^{m} b_i(x)\,\partial_{x_i} + \sum_{j=1}^{n} \beta_j(y)\,\partial_{y_j}.
\end{align*}
Using the coordinate formula for the Lie bracket of vector fields, we obtain
\begin{align*}
[A,B]
=
\sum_{i=1}^{m}\bigl(A(b_i)-B(a_i)\bigr)\,\partial_{x_i}
+
\sum_{j=1}^{n}\bigl(A(\beta_j)-B(\alpha_j)\bigr)\,\partial_{y_j}.
\end{align*}
Because $a_i$ and $b_i$ depend only on the $G$-coordinates and $\alpha_j$ and $\beta_j$ depend only on the $H$-coordinates, the mixed derivatives vanish. Therefore the $G$-component of $[A,B]$ is $[\widetilde X_1,\widetilde X_2]$ and the $H$-component is $[\widetilde Y_1,\widetilde Y_2]$, so
\begin{align*}
[A,B]_{(g,h)}=([\widetilde X_1,\widetilde X_2]_g,[\widetilde Y_1,\widetilde Y_2]_h).
\end{align*}
Evaluating at $(e_G,e_H)$ gives
\begin{align*}
[A,B]_{e_K}=([\widetilde X_1,\widetilde X_2]_{e_G},[\widetilde Y_1,\widetilde Y_2]_{e_H}).
\end{align*}
[/step]