Specific numerical scheme does not converge In my numerical analysis class we are studying numerical solutions to ODEs and we are looking at the equation $y'=-y$ and we want to test the explicit midpoint rule $y_{n+1}=y_{n-1}+2hf(y_n)$ where $f$ is the RHS of $y'=-y$. The initial condition is $y_0=1$ and we take stepsize $h=0.2$ and we want to apply the scheme on the interval $[0,20]$, and the index is interpreted as $y_n \approx y(t = nh)$. And we are asked to take $ y_1 = e^{-h} $. I tried coding the scheme in MATLAB but got the graph of something that oscillates after a certain point and does not look like the true solution at all, but in the smaller interval it does match $e^{-x}$. I have attached the figures below. Can someone please explain why this is happening? I thank all helpers.


 A: On the name of the method
The method is not what is generally known as the explicit midpoint method. For $y'=f(y)$ that would be the one-step method
$$
y_{n+1}=y_n+hf\Bigl(y_n+\tfrac12hf(y_n)\Bigr).
$$
The method described in the question is a multi-step method, while it is also a generalization of the midpoint quadrature method, it is better known as central Euler or two-step Nyström method.
General shape of the leading error terms
Using the definitions and notations of Deriving the central Euler method and intuition, the leading order error terms
can be written as
$$
y_k-y(x_k)=h^2(a_k(x_k)+(-1)^kb(x_k))
$$
where the coefficient functions are solutions of the differential equations
\begin{align}
 a'(x)&=f'(y(x))a(x)-\frac16y'''(x), & 0&=a(x_0)+b(x_0),\\
 b'(x)&=-f'(y(x))b(x), & \frac{y_1-y(x_1)}{h^2}&=a(x_1)-b(x_1).
\end{align}
Applied to the given example
For the test case
$$
y'=f(y)=-y,~~ y(0)=1 \implies y(x)=e^{-x}, ~~y'''=-e^{-x}=-y
$$
this results in
\begin{align}
a'(x)&=-a(x)+\frac16e^{-x}&\implies a(x)&=a(0)e^{-x}+\frac16xe^{-x}\\
b'(x)&=b(x)&\implies b(x)&=b(0)e^{x}.
\end{align}
As per the task $y_1=y(x_1)=e^{-h}$, so that the coefficients $a(0),b(0)$ are obtained by solving the linear system
\begin{align}
0&=a(0)+b(0)\\
0&=a(0)e^{-h}+\frac16he^{-h}-b(0)e^{h}\\[1em]
\hline
a(0)=-b(0)&=-\frac{he^{-h}}{12\cosh(h)}
\end{align}
This means that the oscillating part of the error grows like $\frac{h^3e^{-h}}{12\cosh(h)}e^x$, which for $h=0.2$ has at $x=20$ the numerical value $259603.7064$, which is in the magnitude range of your plot. More precisely,
$$
y_k=e^{-x_k}+\frac{h^2}6x_ke^{-x_k}+\frac{h^3}{12\cosh(h)}(e^{-x_k}-(-1)^{k}e^{+x_k})+O(h^4).
$$
Plotting the error coefficients
The error evolution on the initial segment $x\in [0,4]$ looks like the first plot in the figure below. The colored graphs are the errors $y_k-y(x_k)$ divided by $h^2$ for different sizes of $h$, while the thin lines are the curves $a(x)\pm b(x)$ that the error are expected to lie on or close-by. This they apparently do with only small deviations resulting from higher order error terms.


A completely different, direct but specialized approach
You could also directly solve the linear recursion
$$
y_{n+1}=y_{n-1}-2hy_n
$$
with its characteristic polynomial
$$
0=q^2+2hq-1=(q+h)^2-(1+h^2)\implies q=-h\pm\sqrt{1+h^2}
$$
so that
$$
y_n=(1+c)(\sqrt{1+h^2}-h)^n-c(-1)^n(\sqrt{1+h^2}+h)^n
$$
$y_1=e^{-h}$ then gives
$$
e^{-h}=\sqrt{1+h^2}-h+2c\sqrt{1+h^2}
\\~\\\implies
c=\frac{1-h+\frac12h^2-\frac16h^3+\frac1{24}h^4+O(h^5)+h-1-\frac12h^2+\frac18h^4+O(h^6)}{2\sqrt{1+h^2}}
\\=-\frac{h^3(1-h)}{12}+O(h^5)
$$
The first term goes to zero, the second term for $h=0.2$ and $n=100$ gives the value
$$
\frac{h^3(1-h)}{12}(\sqrt{1+h^2}+h)^{n}\approx 2.2\cdot10^5 
$$
A: I am not sure how much detail you are looking for.  But this is a very good example of when and how numerical methods can fail.  As the value of $x$ increases, errors increase with each step.  For a number of steps the errors remain small, but after a time, the numerical solution deviates and each steps then overshoots by larger and larger errors.  This would happen even if each calculation was perfectly executed.  However, the issue can be compounded by imperfections in computer arithmetic which can increase the unreliability of the numerical method when stretched beyond its limitations.
That is why when handling a numerical method, you need to study the size of the error created at each step and how it changes as you add more steps.
I appreciate this answer is informal, but hopefully of some use.
