I am calculating numerical solutions of time-dependent Schroedinger equation
$\frac{d\Psi}{dt} = - i H \Psi$
where $\Psi$ is an $N$-element complex vector and $H$ is an $N \times N$ complex matrix, which is ``almost normal''. That is, $H = H_0 + i D$, where $H_0$ is normal (and often Hermitian), $D$ is Hermitian but does not commute with $H_0$, and $||D|| \ll ||H_0||$. The dimension $N$ is betwen 30 and 200.
So far, I have been using a 4-th order Runge-Kutta specialised for linear ODEs (taken from Zingg and Chisholm). However, the accuracy is not as good as I'd like, probably due to strong oscillations arising from the $H_0$ part. Is there a better method out there? I'd prefer one which would not require the calculation of $e^{-iH\tau}$ ($\tau$ - time-step of the ODE solver).
EDIT: more data about the problem:
The eigenvalues of $H_0$ are literally random, as $H_0 = h + diag(dE_i)$, where $dE_i ~ N(0,\sigma)$ or $dE_i \sim L(\alpha,\sigma)$. Matrix $h$ has eigenvalues of the order of unity, and $L(\alpha,\sigma)$ is the symmetric Levy alpha-stable distribution (each diagonal value of $H_0$ is perturbed by an independent distr.) with scale parameter $\alpha \in (0,2)$ and strength $\sigma$. $D$ has eigenvalues of the order of 0.1. I choose the time step to be of the order of $0.01/E_\text{max}$, where $E_\text{max}$ is the maximum absolute value of real or imaginary part of any eigenvalue of $H$.
So in short, $H$ can be nasty :(
