The equation and its meaning:
Consider two sets $(A)_{l=0,...,m_a},$,$(B)_{l=0,...,m_b}$ of hermitian matrices and a set of positive semidefinite matrices $(C)_{l=0,...,m_c}$. Each matrix has the dimension $n \times n$. These form the ODE: $$ y^\prime(t) = \sum_{i=0}^{m_a} (A_iy+yA_i^\dagger) + \sum_{i=0}^{m_b}(B_iy-yB_i^\dagger) + \sum_{i=0}^{m_c}C_iyC^\dagger_i. \qquad(*) $$ Note that y(0) is a hermitian $n\times n$ matrix too.
Eqs. like $(*)$ describe the norm-preserving evolution of a linear operator with respect to the anti-commutator (1st), commutator (2nd) and "basis-transform" (3rd). Eqs. like this occur often in Quantum mechanics, see "master equations" for more information.
The problem with vectorization of the matrix matrix form:
The most common ODE solver operates on vectors. So say x(t) is a vector that contains all columns of y(t) then the ODE becomes with the Kronecker product $\otimes$: $$ x^\prime(t) = \left(\sum_{i=0}^{m_a} (1 \otimes A_i + A_i^*\otimes 1) + \sum_{i=0}^{m_b}(1 \otimes B_i - B_i^*\otimes 1) + \sum_{i=0}^{m_c}C_i^*\otimes C_i\right)x \equiv Xx(t). $$ But X is a full $n^2 \times n^2$ matrix! Note that the Jacobi-Matrix coincides with X. So each call to $Xx(t)$ is asymptotically worse than the matrix-matrix multiplication above. Thus, I'm looking for an algorithm to solve the ODE in the matrix form (*) for $y(t_{i=0,...,m_t})$.
A few additional properties
The eigenvalues of the matrices may vary over several magnitudes,
The dimension of $n$ is never bigger than $10^3$, usually it is between $10^1 -10^2$,
The matrices $A_i,B_i,C_i$ may be sparse as single ones, but the sum, i.e. X, is not sparse,
The eigenvalues can be negative,
Stiffness and time-dependent scaling of the matrices, i.e. $A_j \to \alpha(t)_j A_j$, can occur.
Precision should be at least $10^{-5}$, but never more than $10^{-8}$
What I've tried
I know ODEPACK, but they rely on the Jacobian, which is unnecessary big here. RK4 fails, except for very, very small step sizes. So, are they solutions which are easy to adapt for this problem, preferable in Fortran?
Related: