You know by Taylor expansions that $u(t+h)=u(t)+hu'(t+h/2)+\frac{h^3}6u'''(t+h/2)+O(h^5)$. Insert the ODE for the first derivative and compare with the numerical method
\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}
where $\tilde U=U+\frac h2f(t,U)$ so that
\begin{align}
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}
Now let $L_1$ be a bound on $\partial_uf$ around the exact and numerical solution, $M_2$ a bound on $u''$ and $M_3$ a bound on $u'''$. Then inserting these difference expansions into the error gives
\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}
Now assuming that $e=O(h^2)$, which one would have to prove via induction, one can push the $O(he^2)$ terms to the other higher order terms to get
$$
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).
$$
I do not know where your given formula comes from, and what it is intended to prove. As an recursive inequality it shows maximally that $e_n$ is bounded, but no positive error order, while the inequality above leads to
$$e_n\le \frac{e^{L_1(nh)}-1}{L_1}\left[\frac{L_1M_2}8+\frac{M_3}8\right]h^2.$$
