[proofplan]
The argument has two steps. First, assuming we already have a commuting local frame $\hat E_1, \ldots, \hat E_k$ for $E$ (i.e., $[\hat E_i, \hat E_j] = 0$), the composition of their flows provides a local chart whose first $k$ coordinate vector fields coincide with the $\hat E_i$. Second, we build such a commuting frame from involutivity: starting with a chart in which $E|_p$ aligns with $\partial_{x_1}|_p, \ldots, \partial_{x_k}|_p$, we find a local basis of $E$ of the specific form $\partial_{x_i} + \sum_{l>k} a_{il}\partial_{x_l}$, and the structure of this basis forces $[\hat E_i, \hat E_j]$ to lie outside $E$ unless it is zero. Involutivity then compels it to be zero.
[/proofplan]
[step:Commuting flows commute as diffeomorphisms when the fields commute]
We record the following ingredient for the construction. Let $V, W$ be smooth vector fields on an open set with $[V, W] = 0$ on a neighbourhood of $p$. Then by the [commuting flows theorem](/theorems/???), there exists $\varepsilon > 0$ such that the flows $\varphi^V_t, \varphi^W_s$ are defined on a common neighbourhood $N$ of $p$ for $|t|, |s| < \varepsilon$ and satisfy
\begin{align*}
\varphi^V_t \circ \varphi^W_s = \varphi^W_s \circ \varphi^V_t \qquad \text{on } N.
\end{align*}
[/step]
[step:Build a chart from a commuting frame by composing flows]
Assume for now that we have smooth sections $\hat E_1, \ldots, \hat E_k \in \Gamma(E|_U)$ on some open $U \ni p$ that (i) form a pointwise basis of $E_q$ for every $q \in U$, and (ii) satisfy $[\hat E_i, \hat E_j] = 0$ on $U$ for every $1 \le i, j \le k$. We construct the Frobenius chart under this assumption.
Extend $\hat E_1(p), \ldots, \hat E_k(p) \in T_pM$ to a basis $\hat E_1(p), \ldots, \hat E_k(p), v_{k+1}, \ldots, v_n$ of $T_pM$ by adjoining vectors $v_{k+1}, \ldots, v_n \in T_pM$. Choose a smooth chart $(\tilde U, \tilde\varphi)$ at $p$ with $\tilde\varphi(p) = 0$ and with coordinate directions $\partial_{\tilde x_l}|_p = v_l$ for $l = k+1, \ldots, n$ (obtained from any initial chart by a linear change of coordinates on $\mathbb{R}^n$).
Let $\varphi^{\hat E_i}_{t}$ denote the local flow of $\hat E_i$. Define
\begin{align*}
F: (-\varepsilon, \varepsilon)^n &\to M \\
(y_1, \ldots, y_n) &\mapsto \bigl(\varphi^{\hat E_1}_{y_1} \circ \cdots \circ \varphi^{\hat E_k}_{y_k}\bigr)\bigl(\tilde\varphi^{-1}(0, \ldots, 0, y_{k+1}, \ldots, y_n)\bigr),
\end{align*}
for some $\varepsilon > 0$ small enough that all flows are defined and the image stays inside $U \cap \tilde U$. Then $F(0) = p$.
[/step]
[step:Compute the partial derivatives of $F$ using commutativity of the flows]
For $j \in \{1, \ldots, k\}$, we show $\partial_{y_j} F = \hat E_j \circ F$. Because the flows $\varphi^{\hat E_i}_{y_i}$ pairwise commute by $[\hat E_i, \hat E_j] = 0$ and the previous step, we may reorder the composition to place $\varphi^{\hat E_j}_{y_j}$ on the outside:
\begin{align*}
F(y) = \bigl(\varphi^{\hat E_j}_{y_j} \circ \varphi^{\hat E_1}_{y_1} \circ \cdots \widehat{\varphi^{\hat E_j}_{y_j}} \cdots \circ \varphi^{\hat E_k}_{y_k}\bigr)\bigl(\tilde\varphi^{-1}(0, \ldots, 0, y_{k+1}, \ldots, y_n)\bigr),
\end{align*}
where the hat indicates omission. Differentiating in $y_j$ while the other variables are held fixed gives, by the defining property of the flow of $\hat E_j$,
\begin{align*}
\partial_{y_j} F(y) = \hat E_j\bigl(F(y)\bigr).
\end{align*}
For $l \in \{k+1, \ldots, n\}$, only the final term $\tilde\varphi^{-1}(0, \ldots, 0, y_{k+1}, \ldots, y_n)$ depends on $y_l$. Evaluating at $y = 0$,
\begin{align*}
\partial_{y_l} F(0) = (d\tilde\varphi^{-1})_0(e_l) = \partial_{\tilde x_l}\big|_p = v_l.
\end{align*}
[/step]
[step:Apply the inverse function theorem and pull back coordinates to obtain the Frobenius chart]
The differential of $F$ at $0$ has columns $\hat E_1(p), \ldots, \hat E_k(p), v_{k+1}, \ldots, v_n$, a basis of $T_pM$, so $dF_0: \mathbb{R}^n \to T_pM$ is a linear isomorphism. By the [inverse function theorem](/theorems/???), $F$ is a local diffeomorphism onto an open neighbourhood $W$ of $p$. Set $\psi := F^{-1}: W \to (-\varepsilon, \varepsilon)^n$ with coordinates $(x_1, \ldots, x_n) := (y_1, \ldots, y_n) \circ \psi$.
From $\partial_{y_j} F = \hat E_j \circ F$ for $1 \le j \le k$, we get $dF_y(\partial_{y_j}|_y) = \hat E_j(F(y))$, so $\partial_{x_j}|_q = \hat E_j(q)$ for every $q \in W$ and $1 \le j \le k$. In particular,
\begin{align*}
E_q = \operatorname{Span}\langle \hat E_1(q), \ldots, \hat E_k(q)\rangle = \operatorname{Span}\langle \partial_{x_1}|_q, \ldots, \partial_{x_k}|_q \rangle \qquad \text{for every } q \in W.
\end{align*}
This is the desired chart, provided we exhibit a commuting frame.
[/step]
[step:Pick a chart in which $E|_p$ aligns with the first $k$ coordinate directions]
It remains to construct a commuting local basis of $E$. Since $E_p \subset T_pM$ is a $k$-dimensional subspace, choose a linear isomorphism $T_pM \cong \mathbb{R}^n$ sending $E_p$ to $\mathbb{R}^k \times \{0\}$. Pulling back a chart through this linear identification, we obtain a smooth chart $(V, \rho)$ at $p$ with coordinates $(x_1, \ldots, x_n)$, $\rho(p) = 0$, such that
\begin{align*}
E_p = \operatorname{Span}\langle \partial_{x_1}|_p, \ldots, \partial_{x_k}|_p \rangle.
\end{align*}
[/step]
[step:Construct a local basis of $E$ of adapted form $\hat E_i = \partial_{x_i} + \sum_{l>k} a_{il}\partial_{x_l}$]
On $V$, the distribution $E$ is a smooth rank-$k$ subbundle, so it has a smooth local basis $Z_1, \ldots, Z_k \in \Gamma(E|_V)$, each expanded in the coordinate basis as $Z_i = \sum_{m=1}^n Z_{i,m}\, \partial_{x_m}$. At $p$, the matrix $(Z_{i,j}(p))_{1 \le i, j \le k} \in \mathbb{R}^{k \times k}$ has columns giving $Z_i(p)$ projected to $\mathrm{Span}\langle \partial_{x_1}|_p, \ldots, \partial_{x_k}|_p\rangle$. Since $E_p = \mathrm{Span}\langle \partial_{x_j}|_p \rangle_{j \le k}$ and the $Z_i(p)$ form a basis of $E_p$, this matrix is invertible at $p$, hence invertible on a smaller neighbourhood $V' \subset V$ of $p$ by continuity of $\det$.
On $V'$, let $B(q) := (Z_{i,j}(q))_{1 \le i, j \le k}^{-1}$ and define
\begin{align*}
\hat E_i: V' &\to TM|_{V'} \\
q &\mapsto \sum_{j=1}^k B(q)_{ij}\, Z_j(q).
\end{align*}
Each $\hat E_i$ is a smooth section of $E|_{V'}$, and by construction the first $k$ components of $\hat E_i$ in the coordinate basis are the Kronecker symbols $\delta_{ij}$. Thus
\begin{align*}
\hat E_i = \partial_{x_i} + \sum_{l = k+1}^n a_{il}(x)\, \partial_{x_l}
\end{align*}
with smooth coefficients $a_{il} \in C^\infty(V')$. The $\hat E_1, \ldots, \hat E_k$ form a local basis of $E|_{V'}$.
[/step]
[step:Show $[\hat E_i, \hat E_j]$ has no $\partial_{x_l}$-component for $l \le k$]
Compute the bracket in the coordinate basis. Using the general formula for the Lie bracket of coordinate-expanded vector fields from [theorem 1520](/theorems/1520): if $X = \sum_m X_m \partial_{x_m}$ and $Y = \sum_m Y_m \partial_{x_m}$, then
\begin{align*}
[X, Y] = \sum_{m=1}^n \left( \sum_{l=1}^n X_l\, \partial_{x_l} Y_m - Y_l\, \partial_{x_l} X_m \right) \partial_{x_m}.
\end{align*}
Applying this with $X = \hat E_i$, $Y = \hat E_j$, the components $\hat E_{i,m} = \delta_{im}$ for $m \le k$ are constants, hence $\partial_{x_l} \hat E_{i,m} = 0$ for $m \le k$. For $m \le k$,
\begin{align*}
[\hat E_i, \hat E_j]_m = \sum_{l=1}^n \hat E_{i,l}\, \partial_{x_l} \hat E_{j,m} - \hat E_{j,l}\, \partial_{x_l} \hat E_{i,m} = 0 - 0 = 0.
\end{align*}
Therefore
\begin{align*}
[\hat E_i, \hat E_j] = \sum_{l = k+1}^n c_{ij}^{l}(x)\, \partial_{x_l}
\end{align*}
for some smooth functions $c_{ij}^l \in C^\infty(V')$; in particular, $[\hat E_i, \hat E_j]$ has no component along $\partial_{x_1}, \ldots, \partial_{x_k}$.
[/step]
[step:Use involutivity to force $[\hat E_i, \hat E_j] = 0$]
By involutivity of $E$ and the fact that $\hat E_1, \ldots, \hat E_k$ form a local basis of $E|_{V'}$, the bracket $[\hat E_i, \hat E_j] \in \Gamma(E|_{V'})$ expands as
\begin{align*}
[\hat E_i, \hat E_j] = \sum_{s=1}^k \mu_{ij}^s(x)\, \hat E_s = \sum_{s=1}^k \mu_{ij}^s \partial_{x_s} + \sum_{s=1}^k \sum_{l=k+1}^n \mu_{ij}^s\, a_{sl}\, \partial_{x_l}
\end{align*}
for smooth coefficients $\mu_{ij}^s$. Comparing with the previous step, the components along $\partial_{x_s}$ for $s \le k$ must match: $\mu_{ij}^s = 0$ for every $s \in \{1, \ldots, k\}$. Substituting back,
\begin{align*}
[\hat E_i, \hat E_j] = 0 \qquad \text{on } V'.
\end{align*}
The frame $\hat E_1, \ldots, \hat E_k$ on $V'$ satisfies the hypothesis of the chart-building step.
[/step]
[step:Assemble the Frobenius chart]
Apply the chart-building step (with $U := V'$ and the frame $\hat E_1, \ldots, \hat E_k$ constructed in the last two steps) to obtain an open neighbourhood $W \subset V'$ of $p$ and a smooth chart $\psi: W \to (-\varepsilon, \varepsilon)^n$ with coordinates $(x_1, \ldots, x_n)$ (we reuse the symbol $x$; these are the new Frobenius coordinates) such that
\begin{align*}
E_q = \operatorname{Span}\langle \partial_{x_1}|_q, \ldots, \partial_{x_k}|_q\rangle \qquad \text{for every } q \in W.
\end{align*}
The slices $\{(x_1, \ldots, x_n) \in W : x_{k+1} = c_{k+1}, \ldots, x_n = c_n\}$ for each fixed $(c_{k+1}, \ldots, c_n)$ are smooth $k$-dimensional submanifolds of $W$ with tangent space at each point equal to $E_q$, so they are integral manifolds of $E$. Therefore $E$ is integrable, completing the proof.
[/step]