If you solve the equation
$$\frac{dz}{dt}=-i\omega z$$
with initial value $z(0)=1$ you get
$$z = e^{-i\omega t}$$
which traces a circle in the complex plane with angular speed $-\omega$. If you try to the Euler method of approximating the solution, it diverges because
$$z_1 = z_0 + h(-i\omega)z_0 = (1-ih\omega)z_0$$
and
$$|z_1| = |1-ih\omega| = \sqrt{1+(h\omega)^2}$$
for $h,\omega >0$. Likewise $|z_{n+1}|>|z_n|$ for each $n\in\mathbb{N}$.
I'm trying to think of ways to fix this, so that $|z_{n+1}| \leq |z_n|$ and preferably with equality. I tried using the second term of the Taylor series so that the approximation didn't merely make a linear approximation, but the same problem persisted.
It occurs to me that a better method might be to approximate this not with a polynomial but with a circle. I recall from calc 3 the osculating circle, but that was for a curve whose parameterization we already knew, so I'm not sure if I can leverage that. Another thought is that rather than the scheme
$$y_{n+1}=y_n+hy'_n$$
which moves from $y_n$ linearly by adding a multiple of $h$, I might want to rotate using complex multiplication, something like $z_{n+1}=z_nr_ne^{i\omega_nh}$ where $r_n\in\mathbb{R}$ would be the radius of the rotation, larger for when larger curvature is appropriate, and $\omega_n$ controlling the speed and direction--again, of course, $h$ serving as the step-size. I probably would need to do something to deal with the issue of locating the center of the rotation.
I know what I've done is not quite right because I'd have to also account for the center of the rotation, but before going down that road I wanted to know if what I'm attempting is even workable. The big hurdle I can't quite figure out is how I might decide $r_n$ and $\omega_n$ at each stage, using only the derivative information. If I understand it correctly, $\frac{dz}{dt}$ gives only the linear approximation of the direction of the complex number, like with a vector.
You will need implicit methods to preserve the radius. For instance the midpoint method will preserve the radius but not the angular speed. The factor there is $\frac{1-ihω/2}{1+ihω/2}$ so that its absolute value is $1$. Also investigate symplectic methods like Verlet, however there the real and imaginary part are treated differently.