The Frobenius condition which is necessary and sufficient for local integrability is that $Du$ and $D(|\nabla u|^2)$ are parallel. So going back through my answer to this question you see that $u$ must be a function whose gradient descent curves are geodesics.

This also tells you that such functions cannot be dense in $C^2$, answering in the negative your question (1). For example, if you start with a function like $$ f(x,y) = x + y^2 $$ then $Df = (1,2y)$ and $D |\nabla f|^2 = (0, 8y)$ can never be made to be parallel with a small $C^2$ perturbation.

The Frobenius condition which is necessary and sufficient for local integrability is that $Du$ and $D(|\nabla u|^2)$ are parallel. So going back through my answer to this question you see that $u$ must be a function whose gradient descent curves are geodesics.

This also tells you that such functions cannot be dense in $C^2$, answering in the negative your question (1). For example, if you start with a function like $$ f(x,y) = x + y^2 $$ then $Df = (1,2y)$ and $D |\nabla f|^2 = (0, 8y)$ can never be made to be parallel with a small $C^2$ perturbation.

Now, supposing we have a function $u$ such that $D u$ is parallel to $D |\nabla u|^2$. Define $v = |\nabla u|^2$. Then $v$ and $u$ have the same level sets, and thus away from the level sets of $u$ (where $v = 0$) we can write $v = v(u)$ (at least locally).

Letting $F(s)$ be a real valued function such that $F'(s) = \frac{1}{\sqrt{1 + v^2(s)}}$, you see then that the function $F\circ u$ satisfies $$ D(F\circ u) = F'(u) Du = \frac{Du}{\sqrt{1 + v^2(u)}} $$ and so away from the critical points of $u$ this gives you the solution to your equation (*).