Unstable solutions when numerically integrating a system of ODE's with increasing step size Given is the following system of linear ordinary differential equations:
$$ a'(t)= -a(t)+0.1b(t)+0.5c(t)+0.1d(t) $$
$$ b'(t)=0.1a(t)-b(t)+0.2c(t)+0.4d(t)$$
$$ c'(t)=0.5a(t)+0.2b(t)-c(t)+0.3d(t)$$
$$ d'(t)=0.1a(t)+0.4b(t)+0.3c(t)-d(t)$$
with initial conditions
$$ a(0)=1, b(0)=c(0)=d(0)=0 $$
I numerically solved this system with the explicit Euler method at a step size of 0.1 and then increased the step size.
Here is the solution at step size 0.1: 

However, with increasing step size, the solution becomes more and more unstable and eventually "explodes". Here is a graph with step size 1.25:

And here with step size 1.35:

Can somebody explain to me why the increase in step size results in such a behavior?
 A: Say we have a general dynamical system of the form $\dot{x} = Ax$ such that we denote Forward Euler by $x_{n+1} = x_n + Ax_n \Delta x$ for stepsize $\Delta x$. Now we want to know what happens if we perturb the problem a bit; $\tilde{x}_0 = x_0 + \varepsilon$. In the discrete setting we can again write $\tilde{x}_{n+1} = \tilde{x}_n + A\tilde{x}_n\Delta x$. Then if we define the perturbation(error) as $\varepsilon_n = x_n - \tilde{x}_n$ we see that $\varepsilon_{n+1} = (I+A\Delta x)\varepsilon_n$. Of course, for stability we would like $|\lambda | < 1$ for the matrix $(I+A\Delta x)$ which indeed depends on your stepsize $\Delta x$!  
A: You simplify this problem by looking at the solutions of $y' = -\lambda y$ which in theory result to $y={\rm e}^{-\lambda x}$.
The Euler integration assumes the same slope $y'$ for the entire step $h$. So it goes like this
$$\tilde{y}  \leftarrow y+h y' = y - h \lambda y = (1-\lambda h) y  $$
So instead of the next step being $$y \leftarrow y {\rm e}^{-\lambda h}$$
For the above you can notice that $$ 1-\lambda h < {\rm e}^{-\lambda h}$$ always (unless $\lambda h=0$) and there are values that make $1-\lambda h <0$. Such values would cause a flip in the sign of $y$ and you will notice them as unstable oscillations.
This limiting case is when $h>\frac{1}{\lambda}$. In your case you have multiple variables and you can form a matrix/vector product for the ODE slope. Take $\lambda$ as the largest eigenvalue to reach the same conclusion.
