The ODE $${d^2x\over dt^2}=-kx$$ can be converted in the system of linear equations as $$\begin{align} {dx\over dt} & =v\\ {dv\over dt} &= -kx\\ \end{align}$$
Using Euler’s method, given $x_n$ and $y_n$ and for the time step $\Delta t$, the next values can be determined as $$\left[ \begin{matrix} x_{n+1}\\ v_{n+1}\\ \end{matrix}\right] = \left[\begin{matrix} 1&\Delta t\\ -k\Delta t&1 \end{matrix}\right] \left[\begin{matrix} x_n\\ v_n\\ \end{matrix}\right].$$
Now the absolute value of the (possibly complex) eigenvalues should be less than $1$ for this algorithm to be stable. But the eigenvalues turn out to be $1\pm i\sqrt{k}\Delta t$ whose absolute values are strictly greater than $1$ for any nonzero time-step $\Delta t$.
So the algorithm should not work for any value of $\Delta t$, however small. But clearly, this is not the case as my programs do come up with (an approximate) solution though.
So where is the flaw in my reasoning?
The Euler method is indeed not stable for this problem, as the stability analysis states. To see this in the case $k=1$, note that the matrix $$ A = \left[\begin{matrix} 1&\Delta t\\ -\Delta t&1 \end{matrix}\right] $$ satisfies $A^TA = (1 + (\Delta t)^2) I$. That is, $\|(x_{n+1}, v_{n+1})\|^2 = (1 + (\Delta t)^2) \|(x_n,v_n)\|^2$.
Run your approximation long enough and watch the solution spiral out to infinity. Specifically, run it to time $T = (\Delta t)^{-1}$.