Let $U\subset R^n$ be open, $F\colon U\to \mathbb{R}^n$ a $C^\infty$ vector field, and $x(t)$ the solution of $$x’ = F(x)$$ with initial condition $x(0) = y$, which we assume defined at least for $t\in [0,1]$. Define $x_0 = y$ and $$x^\delta_{n+1} = x^\delta_n + \delta F(x^\delta_n),$$ which corresponds to applying Euler’s method with step $\delta$. We know from Euler’s method that (if $\delta$ is small enough) $x^\delta_n - x(n\delta)$ is bounded by $C\delta$ where $C$ is a constant depending only on the bounds of the norms of $F$ and $DF$ in $U$ (the constant is something like $C_2(e^{(n+1)\delta C_1} -1)$ where $C_1, C_2$ are bounds on $F$ and $DF$ respectively). This implies that $x_n^{1/n} \to x(1)$ as $n\to \infty$.
What I can’t find in the literature is the fact that similar estimates hold for higher order derivatives with respect to the initial condition. That is, regarding $x_n^{\delta}$ and $x(t)$ as functions of the initial condition $y$, there are constants $C_k$ and $\delta_k$ depending on bounds on the norms of the derivatives of $F$ up to order $k+1$ in $U$, such that for $\delta<\delta_k$ and $n< 1/\delta$ one has $$\|D^k_y x_n^\delta - D^k_y x(n\delta)\| \leq C_k\delta.$$ The proof seems to be somewhat cumbersome. Does anyone know a reference where this is already done?
