Sabes por las expansiones de Taylor que $u(t+h)=u(t)+hu'(t+h/2)+\frac{h^3}6u'''(t+h/2)+O(h^5)$ . Inserta la EDO para la primera derivada y compárala con el método numérico \begin {align} u(t+h)&=u(t)+hf(t+h/2, u(t+h/2))+ \frac {h^3}6u'''(t+h/2)+O(h^5) \\ U_{+1}&=U+hf(t+h/2, \tilde U) \\ [1em] \hline u(t+h)-U_{+1}&=u(t)-U+h \left [f(t+h/2, u(t+h/2))-f(t+h/2, \tilde U) \right ]+ \frac {h^3}6u'''+... \end {align}
donde $\tilde U=U+\frac h2f(t,U)$ para que \begin {alinear} f(t+h/2, u(t+h/2))-f(t+h/2, \tilde U)&= \partial_uf (t+h/2, u(t+h/2)) \left [u(t+h/2)- \tilde U \right ]+O(|u-U|^2) \\ u(t+h/2)- \tilde U&=u(t)-U+ \frac h2 \left [f(t,u(t))-f(t,U) \right ]+ \frac {h^2}8u''(t) \\ f(t,u(t))-f(t,U)&= \partial_uf (t, u(t)) \left [u(t)-U \right ]+O(|u-U|^2) \end {align} Ahora dejemos que $L_1$ sea un límite para $\partial_uf$ en torno a la solución exacta y numérica, $M_2$ un límite para $u''$ y $M_3$ un límite para $u'''$ . A continuación, insertando estas expansiones de diferencias en el error se obtiene \begin {align} |f(t,u(t))-f(t,U)|& \le L_1e+O(e^2) \\ |u(t+h/2)- \tilde U|& \le e+ \frac h2 \bigl [L_1e+O(e^2) \bigr ]+ \frac {h^2}8M_2 \\ |f(t+h/2, u(t+h/2))-f(t+h/2, \tilde U)|& \le L_1 \Bigl [e+ \frac h2 \bigl [L_1e+O(e^2) \bigr ]+ \frac {h^2}8M_2 \Bigr ]+O(e^2) \\ e_{+1}& \le e+h \Biggl [L_1 \Bigl [e+ \frac h2 \bigl [L_1e+O(e^2) \bigr ]+ \frac {h^2}8M_2 \Bigr ]+O(e^2) \Biggr ]+ \frac {h^3}6M_3+.. \end {align} Suponiendo ahora que $e=O(h^2)$ que habría que demostrar por inducción, se puede empujar la $O(he^2)$ a los otros términos de orden superior para obtener $$ e_{n+1}\le \Bigl[1+L_1h+\frac{(L_1h)^2}2\Bigr]e_n+\frac{h^3}8L_1M_2+\frac{h^3}6M_3+O(h^4). $$
No sé de dónde viene tu fórmula dada, y qué pretende demostrar. Como una desigualdad recursiva muestra como máximo que $e_n$ está acotado, pero sin orden de error positivo, mientras que la desigualdad anterior conduce a $$e_n\le \frac{e^{L_1(nh)}-1}{L_1}\left[\frac{L_1M_2}8+\frac{M_3}8\right]h^2.$$
