Let $\mathbf{u}\colon \mathbb{R}\to \mathbb{R}^{n}$, and suppose $\mathbf{u}$ satisfies $\mathbf{u}' = \mathbf{f}(\mathbf{u}, t)$.
To first order, Euler's method says \begin{align*} \mathbf{u}(t+h) = \mathbf{u}(t) + h \mathbf{f}(\mathbf{u}, t) + \mathcal{O}(h^2) \end{align*} and to second order \begin{align*} \mathbf{u}(t+h) = \mathbf{u}(t) + h \mathbf{f}(\mathbf{u}, t) + \frac{h^2}{2} \left[\nabla_{\mathbf{u}}\mathbf{f}\cdot \mathbf{f} + \frac{\partial \mathbf{f}}{\partial t} \right] \end{align*} However, I have not been able to derive the Taylor expansion to arbitrary order, and even the notation seems to be ambiguous and problematic.
What is the formula for the Taylor expansion of $\mathbf{u}$ to arbitrary order when $\mathbf{u}' = \mathbf{f}$?