Local Truncation Error
You can also go the differentiation route, with $x=x_k$, $x\pm h=x_{k\pm 1}$ and $u_j= u(x_j)$ method step can be related to the the central difference quotient. With the Taylor expansion of $u$ we get
$$
u(x\pm h)=u(x)\pm hu'(x)+\frac{h^2}{2}u''(x)\pm\frac{h^3}6u'''(x)+...
\\~\\
\implies
\frac{u(x+h)-u(x-h)}{2h}=u'(x)+\frac{h^2}{6}u'''(x)+O(h^4)
$$
so that
$$
\frac{u(x+h)-u(x-h)}{2h}-f(x,u(x))=O(h^2)
$$
making this linear multi-step method a second order explicit method.
Global Error Approximation
The distance $d_k=u_k-u(x_k)$ of the numerical method solution $u_k$ to the exact solution $u(x_k)$ evolves as
\begin{align}
u_{k+1}&=u_{k-1}+2hf(u_k)\\
u(x_{k+1})&=u(x_{k-1})+2hf(u(x_k))+\frac{h^2}{6}u'''(x_k)+O(h^4)\\[1em]
\hline
d_{k+1}&=d_{k-1}+2h(f(u_k)-f(u(x_k))-\frac{h^3}3u'''(x_k)+...\\
&=d_{k-1}+2hf'(u(x_k))d_k-\frac{h^3}3u'''(x_k)+...
\end{align}
where the omitted terms are of higher order assuming that $hd_k^2$ is $O(h^4)$ or smaller.
Linear difference equation for the error terms: To get an easily solved problem, first assume that $h$ is small enough so that $u,f'(u),u'''$ are slowly changing over several steps. Then we can set them constant in the above recursion. The now linear recursion
$$
d_{k+1}=d_{k-1}+2hf'd_k-\frac{h^3}3u'''
$$
has for $f'\ne 0$ a solution $d_k=Aq^k+B(-q)^{-k}+C$ where $q>0$ solves
$$
q^2-2hf'q+(hf')^2=1+(hf')^2\implies q=hf'+\sqrt{1+(hf')^2}=e^{hf'+O(h^3)}
$$
The resulting form of the error terms is, up to higher order terms,
$$d_k=Ae^{khf'}+(-1)^kBe^{-khf'}+C.$$
Differential equations for the error components: Translating the terms back into general functions, $d_k=a(x_k)+(-1)^kb(x_k)+c(x_k)$, we identify functions and their differential equations as
\begin{align}
a(x)&\simeq Ae^{(x-x_k)f'(u(x_k))}&\implies a'(x)&=f'(u(x))a(x),\\
b(x)&\simeq Be^{-(x-x_k)f'(u(x_k))}&\implies b'(x)&=-f'(u(x))b(x),\\
c(x_{k-1})-c(x_{k-1})&=2hf'c(x_k)-\frac{h^3}3u'''&\implies c'(x)&=f'(u(x))c(x)-\frac{h^2}{6}u'''(x)
\end{align}
with $u'''(x)=f''(u(x))[f(u),f(u)]+f'(u)^2f(u)$. The initial values are zero for the trend $c$ and account for the error in the first step in the oscillating parts $a,b$.
\begin{align}
c(x_0)&=0, \\
a(x_0)+b(x_0)&=0, \\
a(x_1)-b(x_1)&=u_1-u(x_1)-c(x_1)\approx e(x_0)h^{p+1}+\frac{h^3}{6}u'''(x_0)
\end{align}
This has as consequence that the trend of the error has order $c(x)=O(h^2)$ while the oscillating parts are of order $a(x),b(x)=O(h^3)$ if $p\ge 2$ for the order of the method for the initial step. If the initial step has only order $p=1$, the oscillating parts will reflect this lower order, their contribution will be of the same scale as the error trend curve $c$.
Experimental Error Order Confirmation
That this method is really of order 2 but rather sensitive to the computation of $u_1$ show the following graph depicting the scaled error of a numerical against an exact solution. That the error graphs converge visibly, in the latter paths alternating between an upper and a lower error function, confirms the order, as else the scaled error graphs would be largely different in scale.
The equation used is of the form $F(u, u')=F(p,p')$ with here $F(u,u')=u'+10\sin(u)$ and $p(t)=\cos(t)$, so that $$u'(t)=f(t,u)=10(\sin(\cos(t))-\sin(u(t))-\sin(t).$$ The error for the method applied with step size $h$ is divided by $h^2e^{5t}$, as the equation is slightly stiff and the error rapidly growing.
The initialization is from top down by the exact value $u_1=p(t_1)$, order 3 Heun, by the explicit midpoint method and lastly by the order insufficient Euler method. To decrease the influence of the first error, the step was computed with 5 method steps of step size $h/5$.
See Hairer/Nørsett/Wanner: Solving ODE I: Non-stiff problems, chapter III.9, where Figure 9.2 about this same method has a similar construction.
Best Answer
First, we get local truncation error.
$x(t_{k+1}) = x(t_k) + hf(t_k,x(t_k)) + \tau_k$
$\tau_k = x(t_{k+1}) - x(t_k) - hf(t_k,x(t_k)) = \frac{h^2}{2}x''(\eta)$. Where $\eta \in (t_k,t_{k+1})$.
Then we get the bound,
$$|x(t_{k+1}) - x_{k+1}| \le (1+hL)|x(t_{k}) - x_{k}| + |\tau_k|$$ $$\le (1+hL)|x(t_{k}) - x_{k}| + \frac{h^2}{2}\max_{s \in (0,T)}|x''(s)|$$ $$\le \frac{1}{1-hL}|x(t_{k}) - x_{k}| + \frac{h^2}{2}\max_{s \in (0,T)}|x''(s)|$$
Where the last step is from geometric series. Maybe it would be helpful if you listed the other results you are talking about, and then we can show that they're equivalent.