Motivated by this computation, I was trying to compute this derivative: $$ \left.\frac{d}{dt}\right|_{t=0}X_{\Phi_t^Y(p)} $$ where $X$ is a vector field on a manifold $\mathcal{M}$ and $\Phi_t^Y$ is the flow of the vector field $Y$, i.e $\Phi_t^Y(p)=\gamma(t)$ for $\gamma$ being a curve on $\mathcal{M}$ with $\gamma(0)=p$, $\gamma'(t)=Y_p$.
I am using the definitions in Agricola & Friedrich's Global Analysis so that, if needed, I am allowed to consider that all manifolds are embedded in $\mathbb{R}^n$ for some $n$. Then I can identify a vector $v \in T_p\mathcal{M}$ with a curve $\gamma$ such that $\gamma'(0)=v$ and $\gamma(0)=p$. Also, the definition of differential of a function $f:\mathcal{M}\to \mathcal{N}$ is a linear map $f_{*,p}:T\mathcal{M}\to T\mathcal{N}$ such that: $$ f_{*,p}(v)=\left.\frac{d}{dt}\right|_{t=0}f\circ \gamma(t) $$ since $f\circ \gamma$ is a curve on $\mathcal{N}$.
I have struggled for hours to do the computations with all rigour and finally I am going to post an answer using coordinates, but somehow I feel this should be much more fundamental. Am I missing something? Is there a coordinate-free way to get to the result only from the definitions?