Sorry for the delay, I have been away for a few days. This calculation is horrible and I don't guarantee that this is the best way to do this. I just did it in coordinates and it seems to work.
We want to show that
$$f_*[X,Y]=[f_*(X),f_*(Y)]$$
We will show this by showing that both sides act the same way on functions on $N$ and so must be the same vector field.
Let $h:N\to\mathbb R$, we have by definition
$$f_*([X,Y])(h)=[X,Y](h\circ f)$$
We use coordinates $x^i$ for $M$ and $y^m$ for $N$. We have
\begin{equation}
\begin{aligned}f_*([X,Y])(h)&=[X,Y](h\circ f)(x)\\
&=X^i\frac{\partial}{\partial x^i}\left(Y^j\frac{\partial}{\partial x^j}(h\circ f)\right)(x)-Y^i\frac{\partial}{\partial x^i}\left(X^j\frac{\partial}{\partial x^j}(h\circ f)\right)(x)\\
&=X^i\frac{\partial}{\partial x^i}\left(Y^j\frac{\partial h}{\partial y^m}(f(x))\frac{\partial y^m}{\partial x^i}(x)\right)-Y^i\frac{\partial}{\partial x^i}\left(X^j\frac{\partial h}{\partial y^m}(f(x))\frac{\partial y^m}{\partial x^i}(x)\right)\\
&=X^i\frac{\partial Y^j}{\partial x^i}\frac{\partial h}{\partial y^m}(f(x))\frac{\partial y^m}{\partial x^i}(x)-Y^i\frac{\partial X^j}{\partial x^i}\frac{\partial h}{\partial y^m}(f(x))\frac{\partial y^m}{\partial x^i}(x)
\end{aligned}
\end{equation}
While going from the other end
\begin{equation}
\begin{aligned}
\phantom a[f_*(X),f_*(Y)](h)&=f_*(X)(f_*(Y)(h))-f_*(X)(f_*(Y)(h))\\
&=f_*(X)\left(Y^j\frac{\partial y^m}{\partial x^j}(x)\frac{\partial h}{\partial y^m}(f(x))\right)-f_*(Y)\left(X^j\frac{\partial y^m}{\partial x^j}(x)\frac{\partial h}{\partial y^m}(f(x))\right)\\
&=X^iY^j\frac{\partial y^m}{\partial x^j}(x)\frac{\partial y^n}{\partial x^i}(x)\frac{\partial^2 h}{\partial y^m\partial y^n}(f(x))+X^i\frac{\partial Y^j}{\partial x^m}\frac{\partial y^m}{\partial x^j}(x)\frac{\partial x^m}{\partial y^n}(f(x))\frac{\partial y^n}{\partial x^i}(x)\frac{\partial h}{\partial y^m}(f(x))\\
&\phantom=-Y^iX^j\frac{\partial y^m}{\partial x^j}(x)\frac{\partial y^n}{\partial x^i}(x)\frac{\partial^2 h}{\partial y^m\partial y^n}(f(x))-Y^i\frac{\partial X^j}{\partial x^m}\frac{\partial y^m}{\partial x^j}(x)\frac{\partial x^m}{\partial y^n}(f(x))\frac{\partial y^n}{\partial x^i}(x)\frac{\partial h}{\partial y^m}(f(x))\\
&=X^i\frac{\partial Y^j}{\partial x^i}\frac{\partial y^m}{\partial x^j}(x)\frac{\partial h}{\partial y^m}(f(x))-Y^i\frac{\partial X^j}{\partial x^i}\frac{\partial y^m}{\partial x^j}(x)\frac{\partial h}{\partial y^m}(f(x))\\
\end{aligned}
\end{equation}
We have now shown that both sides are the same. In this derivation I have used the chain rule as well as the Jacobian, some cancellations and the fact that $\frac{\partial y^m}{\partial x^j}(f(x))\frac{\partial x^k}{\partial y^m}(x)=\delta^k_j$.
Suppose $X \in \mathfrak{X}(M), Y \in \mathfrak{X}(N)$ we can define a vector field $X \oplus Y : M \times N \to T(M \times N)$ on product manifold $M \times N$ as
$$
(X \oplus Y)_{(p,q)} = (X_p,Y_q)
$$
under natural identification of $T_{(p,q)}(M \times N)$ with $T_p M \oplus T_qN$ (by isomorphism $\alpha : T_{(p,q)}(M \times N) \to T_pM \oplus T_qN$ defined as $\alpha (v) = (d\pi_M(v), d\pi_N(v))$, one can show that it is a smooth vector field on the product manifold.
So, wlog, given $X \in \mathfrak{X}(M_1)$ it can be checked that for any $X_j \in \mathfrak{X}(M_j)$ for $j=2,\dots,k$, the resulting product $X \oplus X_2 \oplus \cdots \oplus X_k$ is $\pi_1$-related to $X$ by the way the product vector field defined. So vector field on product manifold that $\pi_1$-related to $X$ is not unique. Of course we can choose $X \oplus \mathbf{0}\, \oplus \cdots\oplus \mathbf{0}$ for convenient.
Since you read Lee's, i want to point out that the construction of product vector field above is in fact an exercise in Lee's Introduction to Smooth Manifold (See Problem 8-17 and more general setting in Problem 8-18). However vector fields on the product manifold that $\pi_1$-related to a vector field $X \in \mathfrak{X}(M_1)$ is not necessarily in form of product vector field.
After read this post, i've come to conclusion that
$\mathfrak{X}(M \times N) \supsetneq \mathfrak{X}(M) \oplus \mathfrak{X}(N)$ (as shown in that answer),
Any vector vector field $V$ in product manifold $M \times N$ is in form of $V= X \oplus Y$ for some $X \in \mathfrak{X}(M)$ and $Y \in \mathfrak{X}(N)$ if and only if $V$ and $X$ are $\pi_M$-related and $V$ and $Y$ are $\pi_N$-related.
In more general setting, we know that for any smooth surjective submersion $F : M \to N$ and $X \in \mathfrak{X}(M)$, the pushforward $F_{*}(X)$ is a well-defined smooth vector field on $N$ that is $F$-related to $X$ is and only if $dF_p(X_p) = dF_q(X_q)$ whenever $p$ and $q$ are in the same fiber. So by applying this to the map $\pi_M : M \times N \to M$ and $\pi_N : M \times N \to N$, we have the following criteria :
Any vector vector field $V \in \mathfrak{X}(M \times N)$ is also in $\mathfrak{X}(M) \oplus \mathfrak{X}(N)$ if and only if $d\pi_M(V_{(p,q)})$ constant on each fiber $\{p\} \times N$ and $d\pi_N(V_{(p,q)})$ is constant on each fiber $M \times \{q\}$.
Best Answer
Let $p \in M$. Due to the normal form of smooth submersions, there are smooth chart $(U_p, y_p)$ of $M$, $(V_{f(p)}, x_{f(p)})$ of $N$, containing $p$ and $f(p)$, respectively, such that:
$$ x_{f(p)}\circ f \circ (y_p)^{-1}(z_1, \ldots, z_m) = (z_1, \ldots, z_n)$$
Then, it isn't hard to prove that for $i \in \left \{ 1, \ldots, n\right \}$ $$\mathrm{d}f_{p}\left ( \dfrac{\partial }{\partial y^i}\bigg|_{p} \right ) = \dfrac{\partial }{\partial x^i}\bigg|_{f(p)}$$
I'm not going to write down the subindex $p$ on the components of the charts, to not oversaturate notation. On the other hand, since $X$ is a vector field, it has a representation in $V_{f(p)}$ of the form
$$ X = \sum\limits_{i=1}^{n}g_i^p \dfrac{\partial }{\partial x^i}\bigg|_{f(p)}$$
Afterwards, just define the vector field $Y^p:U_p \to TM$ by
$$Y^p = \sum\limits_{i=1}^{n}\left ( g_i^p \circ f\right ) \dfrac{\partial }{\partial y^i}\bigg|_{p}$$
It is straightforward that $Y^p$ is a smooth vector field and $$\mathrm{d}f_p\left ( Y^p(q)\right ) = X(f(q)) \quad \forall q \in U_p$$
Therefore, we have constructed a family of local vector fields that satisfies the required condition. Now, consider a partition of unity $\left \{ \xi_p\right \}_{p \in M}$ subordinate to $\left \{ U_p\right\}_{p \in M}$ and define $$ Y = \sum\limits_{p \in M}\xi_p Y^p $$
Finally, $Y$ is a global smooth vector field that is $f$-related to $X$.