Criterion for finding invariant sets in continuous dynamical system I'm reading some handouts of a course on dynamical systems, which focuses largely on autonomous systems of ODE's in Euclidean space (i.e. solutions of $\begin{cases} \dot{\mathbf x} = F(\mathbf x) \\ \mathbf x(0) = \mathbf x_0\end{cases}$ where $\mathbf x : \mathbb R \to \mathbb R^n$ and $F : \mathbb R^n \to \mathbb R^n$ is a $C^1$ vector field).
At some point the following statement is given (without proof):

Let $I : \mathbb R^n \to \mathbb R$ be a $C^1$ scalar function and $c \in I(\mathbb R^n)$. Define $\dot I$ as $\dot I(\mathbf x) = \langle \nabla I(\mathbf x), \, F(\mathbf x) \rangle$ (the derivative of $I$ along trajectories). Then, if $\dot I|_{I^{-1}(c)} \equiv 0$, the level set $\{\mathbf x \in \mathbb R^n : I(\mathbf x) = c\}$ is invariant.

Now, it turns out that this is not quite true. Consider, for example, any system with $n = 3$, and let $I(x, \, y, \, z) = x^2 + y^2$, $c = 0$. Then the level set $I^{-1}(c)$ is $\{(0, \, 0, \, z) : z \in \mathbb R\}$ (that is, the $z$-axis). On the other hand, $\dot I(x, \, y, \, z) = 2x\dot x + 2y\dot y$ and its restriction to the $z$-axis is identically $0$, independently of $F$. This would imply (were the statement true) that the $z$-axis is always invariant, which is obviously false.
Or, if you prefer simplicity, just take the one-dimensional case with $I(x) = x^2$ and, again, $c = 0$. $\dot I(x) = 2x\dot x$, which is $0$ for $x = 0$, thus the set $\{0\}$ would have to be invariant (i.e. the point $0$ would have to be fixed), and a trivial counterexample is $F(x) = 1$.
Since this lemma is used extensively throughout the course (for instance, it is used to "prove" the invariance of the $z$-axis in the context of Lorenz equations, by means of the same aforementioned $I(x, \, y, \, z) = x^2 + y^2$), I would like to know if it can be patched somehow, maybe by adding some hypotheses or modifying the statement slightly in order to make it true, so that it still applies to the Lorenz case. I tried a bit but had no luck.
Thanks!

Edit. I thought about it some more and realized that the hypothesis $\nabla I(\mathbf x) \ne 0$ for all $\mathbf x \in I^{-1}(c)$ is necessary. This looks sensible: the condition $\langle I, \, F \rangle = 0$ is saying that $F$ is tangent to the surface $\{I = c\}$ at every point, unless $\nabla I = 0$.
Yet, after a while I think I came up with a counterexample to that as well. Pick $n = 2$ and $F(x, \, y) = (1, \, 1 + y^2)$. Of course $F$ is $C^1$, since $J_F(x, \, y) = \begin{pmatrix} 0 & 0 \\ 0 & 2y \end{pmatrix}$. Now taking $I(x, \, y) = y$ and $c = 0$ we have $\nabla I(x, \, 0) = (0, \, 1) \ne 0$. However, the solutions can be found explicitly and are of the form $\mathbf x(t) = (x_0 + t, \, \tan t)$ for $t \in \left(-\frac{\pi}{2}, \, \frac{\pi}{2}\right)$. So the $x$-axis is not invariant. The only issue might be that we are not in the hypotheses of global existence of the solutions. Nevertheless, intuitively, we are only interested in a neighborhood of the $x$-axis, so we don't really care about what's going on further away. In other words, it should be possible to replace $F$ with some other vector field $\tilde F$ when, say, $|y| > 1$, in such a way that there is a $C^1$ junction in correspondence of the lines $y = \pm 1$ and so that global existence is guaranteed. I believe that, with some patience, one could exhibit an explicit example.

Edit 2. The above counterexample doesn't work, since $\dot I = 1 \ne 0$ (I came up with it thinking the tangent had slope $0$ at $t = 0$...).
In fact, it turned out, with the help of a friend, that the statement becomes true with the additional requirement that $\nabla I \ne 0$. I'll write a proof of it in the answers, for future reference.
 A: Here is a proof that adding the condition $\nabla I(\mathbf x) \ne 0$ if $I(\mathbf x) = c$ makes the statement true. In fact, a stronger statement will be proved:

Let $M \subseteq \mathbb R^n$ be a nonempty, $k$-dimensional differentiable manifold, and suppose that $F(\mathbf x)$ is tangent to $M$ at every point $\mathbf x \in M$. Then, $M$ is invariant.

The lemma stated in the question follows from the latter due to the level set $I^{-1}(c)$ being an $(n - 1)$-dimensional differentiable manifold (this follows from the implicit function theorem when $I$ has nonzero gradient).
Now to the proof. Fix $\mathbf x_0 \in M$ and let $\varphi: U \longrightarrow V$, where $U$ is a neighborhood of $0$ in $\mathbb R^n$ and $V$ is a neighborhood of $\mathbf x_0$, be a diffeomorphism such that $\varphi(0) = \mathbf x_0$ and $\varphi(E \cap U) = M \cap V$, where $E = \mathrm{Span}\{e_1, \, \dots, \, e_k\}$. Any trajectory $t \mapsto \mathbf x(t)$ whose support lies entirely inside $V$ can be mapped to $U$ via $\varphi^{-1}$, i.e. we can define $\mathbf y(t) = \varphi^{-1}(\mathbf x(t))$. Then, $\mathbf y$ satisfies $$\dot{\mathbf y} = J_{\varphi^{-1}}(\mathbf x)\dot{\mathbf x} = (J_{\varphi}(\mathbf y))^{-1}F(\varphi(\mathbf y)) =: G(\mathbf y).$$ Now, whenever $\mathbf y \in E$ (which means $\varphi(\mathbf y) \in M$), $F(\varphi(\mathbf y))$ belongs to the tangent space to $M$ at $\mathbf x 
= \varphi(\textbf y)$, $T_{\mathbf x}M = J_{\varphi}(\mathbf y)(E)$. Thus $F(\mathbf x) = J_{\varphi}(\mathbf y)\mathbf v$ for some $\mathbf v \in E$, and $G(\mathbf y) = (J_{\varphi}(\mathbf y))^{-1}J_{\varphi}(\mathbf y)\mathbf v = \mathbf v$. If we manage to show that the trajctories of $\dot{\mathbf y} = G(\mathbf y)$ starting at any $\tilde{\mathbf y} \in E$ lie (locally) in $E$ (that is, $E$ is invariant for $G$) we are done, since these trajectories are one-to-one with the trajectories of $\dot{\mathbf x} = F(\mathbf x)$ starting at $\tilde{\mathbf x} \in M$ via $\varphi$, and $\varphi$ maps points of $E$ to points of $M$.
To do this, consider any solution of $\dot{\mathbf y} = G(\mathbf y)$ with initial value $\mathbf y(0) = \tilde{\mathbf y}$. If we let $\tilde{\mathbf y} = (\tilde{\mathbf z}, \, 0)$ (where the $0$ is $(n - k)$-dimensional) and, for $\mathbf z \in \mathbb R^k$, $G((\mathbf z, \, 0)) = (H(\mathbf z), \, 0)$ (recall that $G(E \cap U) \subseteq E$), we can consider the system $\begin{cases} \dot{\mathbf z} = H(\mathbf z) \\ \mathbf z(0) = \tilde{\mathbf z} \end{cases}$ in $\mathbb R^k$. Due to the regularity of $H$, local existence and uniqueness of the solution hold. But now, if $\mathbf z(t)$ is the local solution, $\mathbf y(t) = (\mathbf z(t), \, 0)$ is a solution to $\dot{\mathbf y} = F(\mathbf y)$ that lies entirely in $E$, and it is in fact the only one.
