In Hamilton's Ricci Flow (by Chow, Lu, Ni, pp. 29-31, see here) they show that a Riemannian manifold $(M^n,g)$ is locally conformally flat iff the Weyl tensor vanishes (when $n\ge 4$) and iff the Cotton tensor vanishes (when $n=3$).
In the proof they correctly say that it suffices to be able to find a $1$-form $X$ which solves $$ \nabla_i X_j=b_{ij}+X_i X_j-\frac{1}{2}|X|^2 g_{ij}, $$ where $b_{ij}:=\frac{1}{n-2}\left(Ric_{ij}-\frac{1}{2(n-1)}Rg_{ij}\right)$.
Then they appeal to Frobenius theorem and state that it suffices to check that $$ \nabla_k c_{ij}-\nabla_j c_{ik}=-R_{kij}^\ell X_\ell,\ \ \ (*) $$ where $c_{ij}:=b_{ij}+X_iX_j-\frac{1}{2}|X|^2 g_{ij}$, and they proceed by substituting $X$ and $\nabla X$ in order to check $(*)$.
Now to me this makes no sense: what does $(*)$ mean, since we don't know the existence of $X$ yet? How can one substitute $X$ and $\nabla X$ to check a condition before knowing that $X$ exists?
How can this be fixed in order to make sense?
Alternatively, any reference for a good proof is appreciated.
Notice that you are only solving the equation for $X_j$'s in a chart. You have n variables and a priori $n^3$ equations. But the system is symmetric with respect to i,j ... so the number of equations reduce to $n^2$ . It takes a little bit more work to show that the system is also symmetric w.r.t. i and k. So eventually, you'll end up with n equations and n variables. Now these are coefficients of a form so you can first solve the equation point-wise for any point $p$ in a nice coordinate system with $g_{ij}(p) = 0$ and $\partial_k g_{ij}(p) = 0$ and then check to see your solution is still functorial. hope this helps a bit.