A possible proof of Cartan's magic formula $$L_X = i_X \circ d+d \circ i_X$$ is to follow the steps:
- Show that two derivations on $\Omega^{\bullet}(M)$ commuting with $d$ are equal iff they agree on $\Omega^0(M)$.
- Show that $L_X$ is a derivation on $\Omega^{\bullet}(M)$ commuting with $d$.
- Show that $i_X \circ d + d \circ i_X$ is a derivation on $\Omega^{\bullet}(M)$ commuting with $d$.
- Show that $L_X f = Xf = i_Xdf+ d i_Xf$ for all $f \in C^{\infty}(M)=\Omega^0(M)$.
I followed the sketch of proof, and the only point where I am stuck is when I have to prove that $L_X$ and $d$ commute.
More precisely, if I write $$L_Xd\omega =\frac{d}{dt}_{|t=0} \phi_t^* d\omega = \frac{d}{dt}_{|t=0} d\phi_t^* \omega = \lim\limits_{t \to 0} d \left( \frac{1}{t} \left( \phi_t^* \omega- \omega \right) \right),$$ is there an argument to permute $d$ and $\lim\limits_{t \to 0}$?
Putting in the limit made it more mysterious. You're just using the fact that for smooth functions mixed partials are equal. That is, $$\frac{\partial}{\partial t}\frac{\partial}{\partial x^j}=\frac{\partial}{\partial x^j}\frac{\partial}{\partial t}.$$
EDIT: In particular, if you write it out in local coordinates, it's computed by taking partial derivatives with the respect to the variables $x^1,\dots,x^n$. Writing $\phi_t^*\omega=\sum f_I(x,t)dx^I$ (where $I=(i_1,\dots,i_k)$ ranges over length-$k$ multiindices), then $d\phi_t^*\omega=\sum \frac{\partial f_I(x,t)}{\partial x^j}\,dx^j\wedge dx^I$, etc.