I'm trying to solve this ODE system:
$$u''(x)=-9000u(x)-5000(\sin(u(x)) + u(x)),\ \ u(0)=\pi,\ u'(0)=0.$$ With $x \in (0,10)$.
I have to solve it with Euler (explicit). So $y_n+1=y_n + kf(t_n,y_n)$. The problem is the A-stability. If I take 200 time-steps, say $ts=200$, then $k=10/200$ and the numerical solution diverge. Then, inspired by the restriction on $k$ (see here) I took $k$ really small, for example $k=1e-8$ and I got the correct numerical solution.
Is there a better formal way to take this $k$?
Your system is conservative with potential function $$V(u)=7000u^2+5000(1-\cos(u)).$$ As the initial point is at rest, it also marks the highest point of the potential, thus the exact solution is constrained to $u(x)\in[\pi,\pi]$ and is a periodic oscillation. Estimating the average influence of the sine term over that interval as something like $0.2u$ gives an approximation by the linear oscillator $$u''+15000u=0$$ with frequency $\omega=100\sqrt{1.5}\approx 122.5$. By the sampling theorem you need at least 2, for any sensible purposes at least 4 samples per period to capture the oscillation somewhat faithfully. This tells you that the sampling step, independent of the method, should be smaller than $\frac{2\pi}{490}\approx 0.012$ or at least $800$ time steps.
Using a balanced first order system close to $$ u'=\ \ ωv\\v'=-ωu $$ you get imaginary eigenvalues that never fall into the stability region of the forward Euler method. Falling back to a more heuristic approach, take the Lipschitz constant $L$ which is close to $ω$, and apply the guideline $Lh\ll0.5$ which suggests to try $h\ll 0.004$, that is more than $2500$ steps for the given interval to get something close to stability. This just means that the numerical approximation does not explode directly from the start in a qualitatively wrong manner.
However, the Euler method spirals outwards in the phase picture, as the level sets of the energy function are more or less concentric circles. This happens with a factor to the radius in the phase plane of about $\sqrt{1+(ωh)^2}$ per step or $\exp(5ω^2h)$ over the full integration interval. To get the radius increase lower than a factor $10$ one needs $5ω^2h<3$ or $h<0.00004$, which are $250000$ time steps.