[proofplan]
We work in a local chart near an arbitrary point and write $X, Y$ in terms of the coordinate basis. Direct computation shows that the second-order partial derivatives that appear in $X(Y(f))$ and $Y(X(f))$ cancel by symmetry of mixed partials, leaving a first-order expression whose coefficients are smooth functions of the coordinates. This exhibits $[X,Y]$ locally as a derivation of $C^\infty$, hence a smooth vector field.
[/proofplan]
[step:Set up a coordinate chart and write $X, Y$ in the coordinate basis]
Fix $p \in M$ and choose a smooth chart $(U, \varphi)$ at $p$ with coordinates $(x_1, \ldots, x_n)$ where $\varphi: U \to \varphi(U) \subset \mathbb{R}^n$ is a diffeomorphism onto an open set. On $U$, the coordinate vector fields $\partial_{x_1}, \ldots, \partial_{x_n}$ form a smooth local frame for $TM|_U$. Write
\begin{align*}
X|_U &= \sum_{i=1}^n X_i\, \partial_{x_i}, & Y|_U &= \sum_{j=1}^n Y_j\, \partial_{x_j},
\end{align*}
with component functions $X_i, Y_j \in C^\infty(U)$. These exist and are smooth because $X, Y$ are smooth sections of $TM$.
[/step]
[step:Expand $X(Y(f))$ and $Y(X(f))$ by the Leibniz rule]
Let $f \in C^\infty(U)$. Applying $Y$ to $f$ gives $Y(f) = \sum_j Y_j\, \partial_{x_j} f \in C^\infty(U)$. Applying $X$ and using the Leibniz rule at each step,
\begin{align*}
X(Y(f)) &= \sum_{i=1}^n X_i\, \partial_{x_i}\!\left( \sum_{j=1}^n Y_j\, \partial_{x_j} f \right) \\
&= \sum_{i=1}^n \sum_{j=1}^n X_i\, (\partial_{x_i} Y_j)\, (\partial_{x_j} f) + \sum_{i=1}^n \sum_{j=1}^n X_i\, Y_j\, \partial_{x_i}\partial_{x_j} f.
\end{align*}
By the symmetric computation with $X$ and $Y$ exchanged,
\begin{align*}
Y(X(f)) &= \sum_{i=1}^n \sum_{j=1}^n Y_j\, (\partial_{x_j} X_i)\, (\partial_{x_i} f) + \sum_{i=1}^n \sum_{j=1}^n Y_j\, X_i\, \partial_{x_j}\partial_{x_i} f.
\end{align*}
[/step]
[step:Cancel the second-order terms by symmetry of mixed partials]
Subtracting, the second-order sums are
\begin{align*}
\sum_{i,j} X_i Y_j\, \partial_{x_i}\partial_{x_j} f - \sum_{i,j} Y_j X_i\, \partial_{x_j}\partial_{x_i} f.
\end{align*}
Since $f \in C^\infty(U)$, [Clairaut's theorem on equality of mixed partials](/theorems/???) gives $\partial_{x_i}\partial_{x_j} f = \partial_{x_j}\partial_{x_i} f$. The component functions $X_i, Y_j$ commute as elements of $C^\infty(U)$, so $X_i Y_j = Y_j X_i$ pointwise. The two sums are therefore equal and cancel:
\begin{align*}
\sum_{i,j} X_i Y_j\, \partial_{x_i}\partial_{x_j} f - \sum_{i,j} Y_j X_i\, \partial_{x_j}\partial_{x_i} f = 0.
\end{align*}
[/step]
[step:Identify the remaining first-order expression as a derivation with smooth coefficients]
After cancellation,
\begin{align*}
[X, Y](f)\big|_U = X(Y(f)) - Y(X(f)) = \sum_{i=1}^n \left[ \sum_{j=1}^n \bigl( X_j\, \partial_{x_j} Y_i - Y_j\, \partial_{x_j} X_i \bigr) \right] \partial_{x_i} f.
\end{align*}
Define
\begin{align*}
A_i: U &\to \mathbb{R} \\
q &\mapsto \sum_{j=1}^n \bigl( X_j(q)\, \partial_{x_j} Y_i(q) - Y_j(q)\, \partial_{x_j} X_i(q) \bigr),
\end{align*}
so that
\begin{align*}
[X, Y](f)\big|_U = \sum_{i=1}^n A_i\, \partial_{x_i} f.
\end{align*}
Each $A_i$ is a finite sum of products of $C^\infty(U)$ functions ($X_j, Y_j$ and their first-order partial derivatives), hence $A_i \in C^\infty(U)$.
[/step]
[step:Conclude that $[X,Y]$ is a smooth vector field on $M$]
For every $q \in U$ and $f \in C^\infty(M)$, the operator
\begin{align*}
[X,Y]_q: C^\infty(M) &\to \mathbb{R} \\
f &\mapsto \sum_{i=1}^n A_i(q)\, \partial_{x_i} f(q)
\end{align*}
is $\mathbb{R}$-linear and satisfies the Leibniz rule at $q$ (being a linear combination of the derivations $\partial_{x_i}|_q$), so $[X,Y]_q \in T_q M$. The coefficient functions $A_i$ are smooth on $U$, so the assignment $q \mapsto [X,Y]_q$ is a smooth section of $TM$ over $U$.
Since $p \in M$ was arbitrary and the definition of $[X,Y]$ is chart-independent, the section $[X,Y]$ is smooth on all of $M$. Therefore $[X,Y] \in \Gamma(TM)$, completing the proof.
[/step]
