Hi,
Two points:
- If Omega_n is an operator, it cannot be given as coefficients. The coefficient should be a function that returns a float/complex number. In this case, there is no need to distinguish between coefficients and Hamiltonian. You can write a Python function that returns HI and give it to the solver.
```
def HI(t, args):
# type the equation of HI here.
mesolve(HI, ...)
```
- For exponentiation of operators, np.exp is wrong. One should use `(-1.j * Omega * times).expm()`, which is same as scipy.linalg.expm
Best, Boxi