In what follows I am assuming (but I'm pretty sure) that your integrable sub-distribution $\mathscr{L}$ is supposed to be the same distribution as $\mathscr{F}$ that you mention later in your question but please correct me if I've misunderstood. Additionally, I'll assume we're working locally, i.e. on sufficiently small open subsets of a fixed point $p$ in $\mathbb{R}^{2m+1}$ (in order to use Frobenius theorem). Oh, and I'll be using the Einstein summation convention, too.
You are right that the Frobenius theorem is key here. Indeed, let $z, y^1,\ldots, y^m$ be the first integrals of $\mathscr{L}$ (i.e. $X(z)=X(y^i)=0$ for all $1\leq i\leq m$ and for all $X\in\mathscr{L}$). These help define a set of coordinates on $\mathbb{R}^{2m+1}$ on $(z, y^1,\ldots,y^m,x^1,\ldots, x^m)$ so that
$$
\mathscr{L}=\text{span}\{\partial_{x^1},\ldots,\partial_{x^m}\},
$$
as you had already noticed. Now, since the first derived system, call it $\mathscr{D}'$, (informally defined as $\mathscr{D}'=\mathscr{D}+\left[\mathscr{D},\mathscr{D}\right]$), then the growth vector of $(m+1,n)$ does indeed guarantee that $\mathscr{D}'=T\mathbb{R}^n$. Since $\mathscr{L}$ has rank $m$, this means that
$$
\mathscr{D}=\mathscr{L}\oplus\text{span}\{Z\}
$$
for some vector field $Z$ with no components from the directions in $\mathscr{L}$. In particular, since $\text{rank}\,\mathscr{D}'=2m+1$ and $\mathscr{L}$ is integrable, then the vector fields defined by $Y_i=[\partial_{x^i},Z]$ for all $1\leq i\leq m$ must not be in the distribution $\mathscr{D}$ (but they are of course in $\mathscr{D}'$). Writting $Z=A\partial_{z}+B^i\partial_{y^i}$ then one can make the argument that the first integrals of $\mathscr{L}$ may be written in such a way so that
\begin{equation*}
\begin{aligned}
Y_1&=\frac{\partial A}{\partial x^1}\partial_z+\frac{\partial B^i}{\partial x^1}\partial_{y^i}=\partial_{y^1}\\
\,&\vdots\\
Y_m&=\frac{\partial A}{\partial x^m}\partial_z+\frac{\partial B^i}{\partial x^m}\partial_{y^i}=\partial_{y^m}.
\end{aligned}
\end{equation*}
This is a system of linear 1st order PDE and we can quickly see a solution to be $A=A(y^1,\ldots,y^m)$ (i.e. no dependence on the $x^i$'s) and $B^i=x^i\tilde{B}^i(y^1,\ldots,y^m)$ for all $1\leq i\leq m$. Thus the distribution $\mathscr{D}$ can be written as
\begin{equation*}
\mathscr{D}=\text{span}\{A\partial_z+x^1\tilde{B}^1\partial_{y^1}+\cdots+x^m\tilde{B}^m\partial_{y^m},\partial_{x^1},\ldots,\partial_{x^m}\}
\end{equation*}
Edit: If I recall correctly, this is essentially a special case of Robert Bryant's doctoral thesis, so you may also find that as an enlightening resource.
