ODE for the differential of a flow on a manifold This question is a consequence of my horrible knowledge in differential geometry. It can be stated as follows.
Consider the solution $y_p(t)$ to the ODE $$\partial_t y_p(t) = X(y_p(t)), \qquad y_p(0) = p$$ where we are given a manifold $M \ni p$ and a vector field $X$ (everything is smooth). This induces a flow $\varphi(t,p) = y_p(t).$ In Euclidian space I can write the equation for the derivative of this flow:
$$
\partial_t D_x\varphi(t,x) = D_x X (\varphi(t,x)) \cdot D_x\varphi(t,x), \qquad D_x \varphi(0,x) = Id.
$$
My question is: how do I write this equation in and "invariant" way on a manifold? 
Crucially, what is the correct generalisation for the matrix product $ D_x X (\varphi(t,x)) \cdot D_x\varphi(t,x)$?
In my not-knowledge of differential geometry I learned that differentiating vector fields it one of the slightly complicated things...
Similar questions have been asked. Here the same question appeared, but the answer regards the determinant of the Jacobian and not the Jacobian, I believe (https://mathoverflow.net/questions/284718/derivative-of-the-flow-for-odes-on-manifolds).
EDIT: In the answer I posted below it appears that given a connection it is possible to write down an ODE for the flow in the classical sense. The following question is still open: 
Is it possible to write an ODE for the differential in a generalized way (without assuming the existence of a connection), such that given a connection the equation reduces to the classical one?
 A: Here are some calculations, which I post as an answer, so that no confusion arises.
Let us suppose we have a Riemannian manifold admits a symmetric connection. Fix $p \in M, q \in T_p M.$ We consider the differential $D_p \varphi_t (q) \in T_{\varphi_t(p)} M$. We see that we have built a vector field along the curve $\varphi_t(p)$. It thus makes sense to differentiate along the curve to find an equation. We have:
$$D_t (D_p \varphi_t(q)) = q^j \frac{\partial}{\partial x^j} [ \partial_t \bar{x}^k(\varphi_t(p)) ] \frac{\partial}{\partial \bar{x}^k} + X^i(\varphi_t(p)) q^j \frac{\partial \bar{x}^\ell (\varphi_t(p))}{\partial x^j} \Gamma^k_{i \ell} \frac{\partial}{\partial \bar{x}^k}.$$
Under the assumption of symmetry this can be rewritten as:
$$
D_t (D_p \varphi_t(p)) = \nabla_{D_{p} \varphi_t(q)} X = \nabla X \cdot D_{p} \varphi_t(q)
$$
where the $\cdot$ indicates a contraction along the correct indeces.
On the other hand if we do not assume symmetry it appears that we find the equation:
$$
D_t (D_p \varphi_t(p)) =  \nabla X \cdot D_{p} \varphi_t(q) + \tau(X, D_{p} \varphi_t(q))
$$
where $\tau$ is the torsion tensor.
