Suppose that you have a system of ODEs $$\dot{x} = F(x), \, x \in \mathbb{R}^n,$$ and $x_0$ is a an equilibrium point of this system, i.e. $F(x_0) = 0$. Suppose you have a locally invariant smooth submanifold $\mathcal{M}$ that passes through $x_0$. By locally invariant I mean that for any point $p \in \mathcal{M}$ there exists an $\varepsilon > 0$ such that $\varphi^{t}(p) \in \mathcal{M}$ for $\lvert t \rvert < \varepsilon$, where $\varphi^t$ is the flow defined by a system of ODEs. I believe that if these conditions are met, then $DF(x_0) (T_{x_0}\mathcal{M}) \subset T_{x_0}\mathcal{M}$; $DF(x_0)$ is a Jacobi matrix for $F(x)$ at point $x_0$.
It feels like I had encountered similar statement somewhere before, but I can't remember the exact reference. I can't find a counter-example or a proof for this statement from the get-go, though I definitely want to crack that puzzle if there is no reference to it. If this statement was true, it would much simplify one of my proofs, so that's the reason why I'm so interested in it.
So, mostly I'm interested in the reference where this (or very similar) statement had been proved. If the statement is false and you know a counter-example, that would be very useful (although a bit heartbreaking).
ADDED LATER
I have a sketch of the proof for this statement (kind of). If we pick any vector $v \in T_{x_0}\mathcal{M}$ there exists a curve $\gamma_v(s)$ such that $\gamma_v'(0) = v, \gamma_v(0) = x_0,$ and which lies in $\mathcal{M}$. If we choose some $\tau$ and apply $\varphi^{\tau}$ to $\gamma_v(s)$ we will obtain some curve $\tilde{\gamma}_v(s)$ that still lies on $\mathcal{M}$ and $\tilde{\gamma}_v(0) = x_0$. Thus tangent vector to this curve has to lie in $T_{x_0}\mathcal{M}$. But this tangent vector is nothing but $DF(x_0)(v)$ and this implies that $\forall v \in T_{x_0}\mathcal{M}$ follows $DF(x_0)(v) \in T_{x_0}\mathcal{M}$.
Since $F(x_0) = 0$, the coordinate derivative $\def\M{\mathcal{M}}D F(x_0)=\partial_i F^j(x_0)\,dx^i\otimes\partial_j$ is a well-defined endomorphism of $T_{x_0} \M,$ defined independently of the coordinates we use to compute it. (One way to see this is via the Riemannian formula $\nabla_i F^j = \partial_i F^j + \Gamma_{ik}^j F^k= \partial_iF^j.$) This means we can perform arbitrary smooth coordinate changes without changing the meaning of $DF$.
Now that we know we're allowed to, choosing the right coordinates makes this problem very easy. We will use submanifold slice coordinates for $\M;$ i.e. local coordinates centred on $x_0$ so that $\M$ coincides locally with the $m-$plane $$\{ x \in \mathbb R^n : x^{m+1} = 0, \ldots, x^n = 0\}.$$
The fact that $\M$ is invariant tells us that $F|_\M$ lies in the same coordinate plane, i.e. $F^{m+1} = 0,\ldots,F^n = 0$ along $\M.$ Thus we clearly have $DF^{m+1}(v) = 0, \ldots, DF^n(v) = 0$ for any direction $v \in T_{x_0} \M$, which tells us exactly that $DF(T_{x_0} \M) \subset T_{x_0} \M.$
Here's an alternate proof making your idea from the OP rigorous. For any $v \in T_{x_0} \M,$ choose a curve $\gamma : (-\epsilon,\epsilon) \to \M$ with $\dot \gamma(0) = v$ and consider the two-parameter map $\sigma(s,\tau) = \varphi^\tau \gamma(s).$ Since $\varphi$ is the flow of $F,$ differentiating this with respect to $\tau$ yields $\partial_\tau \sigma(s,0)=F(\gamma(s))$, which we can differentiate with respect to $s$ to find $$\partial_s \partial_\tau \sigma(0,0) = DF(\dot \gamma(0))=DF(v).$$ On the other hand, since $F(x_0) = 0$ and $\M$ is invariant, we know that for any $\tau$ the curve $\varphi^\tau \gamma$ starts at $x_0$ and stays in $\M$; so $\partial_s \sigma(0,\tau) \in T_{x_0} \M$ for each $\tau$ and thus $\partial_\tau \partial_s \sigma(0,0) \in T_{x_0} \M.$ Partial derivatives commute, so we are done.