Maintain unitary time evolution for a nonlinear ODE

Crank-Nicholson is usually used for linear ODE but can also be extended to a nonlinear ODE

$$dot z = f(z)$$

like so:

$$frac{z_{n + 1} - z_n}{delta t} = frac{f(z_{n + 1}) + f(z_n)}{2}.$$

Another similar alternative is the implicit midpoint rule:

$$frac{z_{n + 1} - z_n}{delta t} = fleft(frac{z_{n + 1} + z_n}{2}right).$$

The Crank-Nicholson and implicit midpoint rules are the same for linear ODE, but for nonlinear ODE, they're different. The reason I mention the implicit midpoint rule is that it's a symplectic integration scheme, whereas the Crank-Nicholson method is not. If you haven't run across this concept before, I highly recommend reading Geometric Numerical Integration by Hairer, Lubich, and Wanner*.

One blunt approach you could take is to approximate $M$ as piecewise constant in each timestep and use the Pade approximant

$$e^{iz} approx frac{1 + frac{1}{2}iz}{1 - frac{1}{2}iz}$$

for $e^{-idelta tM}$. This approximation is especially nice because it wraps the real line into the complex unit circle. When you then apply it to matrices, it takes things with real eigenvalues (like Hermitian matrices) and maps them into unitary matrices. You then get the following scheme:

$$begin{align} A_{n + 1} & = A_n cdot left(I - frac{1}{2}idelta t M_nright)left(I + frac{1}{2}idelta tM_nright)^{-1} M_{n + 1} & = A_{n + 1}^dagger H(t_{n + 1})A_{n + 1} end{align}$$

If $A_0$ is a unitary matrix, then every successive iterate will be a unitary matrix (at least in real arithmetic). It's only a first-order accurate integration scheme, but the important part is that it preserves the structure of your system. Hopefully this is a good starting point from which you can derive more accurate schemes if need be.

*As Lutz Lehmann pointed out in the comments, your problem isn't a Hamiltonian system as such, so what I said about symplectic integrators is a bit of a non-sequitur. I mentioned this more just for general knowledge. Your particular problem has non-trivial geometric structure and so do Hamiltonian systems. If you often find yourself solving these kinds of problems then the knowledge you acquire about methods for Hamiltonian systems will transfer well to other types.