Solution (source code)

= Solution

Put $e_n=y(x_n)-y_n$. Subtracting the numerical update from the exact update with its <local truncation error> gives
$$
e_{n+1}=e_n+h\theta[f(x_{n+1},y(x_{n+1}))-f(x_{n+1},y_{n+1})]+h(1-\theta)[f(x_n,y(x_n))-f(x_n,y_n)]+hT_n.
$$
The <Lipschitz condition> in the solution variable and the <triangle inequality> imply
$$
\boxed{(1-L\theta h)|e_{n+1}|\leq[1+L(1-\theta)h]|e_n|+h|T_n|.}
$$
The assumed positive denominator permits division. With $K=[1+L(1-\theta)h]/(1-L\theta h)$ and $e_0=0$, repeated substitution of this <Lipschitz one-step error bound> yields
$$
|e_n|\leq\frac h{1-L\theta h}\sum_{j=0}^{n-1}K^{n-1-j}|T_j|
\leq\boxed{\frac h{1-L\theta h}\frac{K^n-1}{K-1}\max_{0\leq j<n}|T_j|.}
$$
Here $K>1$ because $L,h>0$, so the displayed <geometric series> quotient is well defined. For implicit steps, the same condition makes the scalar update map a <contraction mapping>, ensuring a unique numerical next value.