Consider the ODE $$y'(t) = \lambda y(t), \quad y(t_0) =y_0.$$ Euler's method $y_{i+1}=y_i+h\lambda y_i $ is stable (meaning that the solution decays or stays constant as $ i \to \infty$) provided that $ |1 + h\lambda| \leq 1 $. This idea can be extended to systems of $n$ dimensions by placing a similar condition on the maximum coefficient instead of just $\lambda$.
This question is about non-linear, $n$ dimensional systems. Under what conditions is Euler's method stable for the system
$$\mathbf{y}'(t)=\mathbf{F}(\mathbf{y}(t)), \quad \mathbf{y}(t_0)= \mathbf{y}_0$$?
I suspect this will involve the Jacobian of the non-linear term. The only resource I've found is on slide 18 of these lecture slides, but the notation isn't explained (what is $l_{k+1}$?). It was probably introduced in a previous lecture. The integral of the Jacobian that appears here is interesting, it looks like some kind of weighted average between the numerical and the true solution. Where does this come from?
Given the linear problem $\vec{y}'=A\vec{y} + \vec{b}$, where $A$ is a diagonalizable n-dimensional matrix, you can find matrices $V,D$ s.t $AV=VD$, with $D=diag${$\lambda_1,...,\lambda_n$}.
Then $\vec{y}'=A\vec{y} + \vec{b} $ iff $ \vec{z}'=D\vec{z} + V^{-1}\vec{b}$. Since we're using Euler's method, the restriction on the timestep is given by the smallest eignevalue of $A$.
If we consider now $\mathbf{y'}=\mathbf{F(y(t))}$,we can linearize $\mathbf{y'(t)})=\mathbf{F(y(t))} + J(\mathbf{y_n})(\mathbf{y}(t) - \mathbf{y_n})$, where $J$ is the jacobian matrix. If $J$ is diagonalizable, then we proceed as done in the linear case (if $J$ is not diagonalizable, then we use Jordan decomposition.
Problems arise because the eigenvalues of $J$ aren't usually enough to describe the behaviour of the exact solution.
We can have several cases: for example $J_f$ could has $\Re(\lambda_i)<0$, but the solution doesn't doesn't decay monotous for every t.
That's what happen if we consider the system $y'=Ay$, with $A=[ ( -\frac{1}{2t} , \frac{2}{t^3} ),( -\frac{t}{2} , -\frac{1}{2t} ) ] $
For more references, see Lambert J. (1991) Numerical methods for ODE systems