It is a known fact that a Manifold $M$ is locally flat if and only if the Riemann curvature tensor (Levi-Civita connection) vanishes. To show the if part is easy, but I am trying to follow the details of the proof given here but I am at a loss in reaching a full proof. I am interested only in a proof of this statement using the Frobenius theorem.
Here is what I have tried: We start with a coframe $\Phi$, it obeys the structure equations ($R=0,~T=0$)
$$ d\Phi + \omega \wedge \Phi =0, \qquad d\omega + \omega \wedge \omega=0 $$ from the second equation it can be seen (Frobenius thm) that the first equation is integrable. Now, since $\Phi$ is a coframe I can ask whether a local rotation matrix $R(x) \in SO(n)$ will be capable of rotating the coframe $\Phi$ to a closed coframe $\Theta$, $d\Theta=0$ (If we manage to do that then Poincaré lemma implies that $\Theta = d \Psi$ locally, and such $\Psi$ are the sought euclidean coordinate functions, we are done). Calculation gives:
$$ d(R \Phi) = d \Theta = (dR - R \omega)\wedge \Phi $$ and this will vanish if the equation $R^{-1}dR = \omega $ is solved by some $R$. I can see that the structural equation implies the integrability of this equation, so a solution $R$ exists, but how can we be sure that such solution will be a matrix in $SO(n)$?
Here is my question: How do we know that there is an $R \in SO(n)$ that solves $dR = R \omega$?
Following Élie Cartan, you want to think of your differential system $dR - R\omega = 0$ on $M\times SO(n)$. This $\mathfrak{so}(n)$-valued $1$-form is integrable, as you said, because of vanishing curvature. The integral manifolds of this differential system will locally be the graphs of functions $R\colon U\to SO(n)$.