Super Adiabatic Evolution numerically I'm doing some adiabatic and super adiabatic evolution using python but without success, the super adiabatic to be more precise. The problem is the way I've manage to write the super adiabatic Hamiltonian, it can be written in function of the eigenvalues and eigenvectors of the Adiabatic Hamiltonian as:
$$ H_{SAD} = \sum_n \left [ E_n|n(t)\rangle\langle n(t)|+\mathrm{i}\hbar\left (|\dot n (t)\rangle \langle n(t)| - \langle n(t)|\dot n (t)\rangle|n(t)\rangle\langle n(t)| \right )\right ]. $$
My two different approaches were a numerical and a symbolic. In the numerical I've created a list of Adiabatic Hamiltonians (so i get one array with this dimensions [ row of Hamiltonian, column of Hamiltonian, time]). The problem in this approach is that I'm not sure how to solve the dynamics of this system that are ruled by the Master Equation:
$$ \dot \rho_{\mathrm{tot}} = -\frac{\mathrm{i}}{\hbar}[H_{\mathrm{tot}}, \rho_{\mathrm{tot}}(t)]. $$
I've found this Python library that solves matrix ODE (odeintw), but a matrix in the shape of mine are not allowed in it, I would need a explicit time dependence to do so:
I've tried to use mesolve from QuTip but the problem is the same, the solver doesn't support a time dependence express as a list of matrices.
This is how I construct my Super Adiabatic Hamiltonian:
import numpy as np
from qutip import *
import matplotlib.pyplot as plt
from odeintw import odeintw

x,y,z, I = sigmax(), sigmay(), sigmaz(), qeye(2)
l = 3
psi_i = tensor(bell_state('11'),ket('1'))
rho0 = ket2dm(psi_i)

Hi = tensor(x,x,I) + tensor(y,y,I)
#Normalized t, s = t/tau
s = np.linspace(0, 1, 100) 

#Adiabatic Hamiltonian
Had = np.zeros([2**l,2**l,len(s)],dtype=complex)

#Super Adiabatic Hamiltonian
Hsad = np.zeros([2**l,2**l,len(s)],dtype=complex)

#Eigenvectors of Had
n = np.zeros([2**l,2**l,len(s)],dtype=complex) 


#Time derivative of n
dn = np.zeros([2**l,2**l,len(s)],dtype=complex)

#Eigenvalues of Had
E = np.zeros([len(s),2**l],dtype=complex)

i = 0 
for t in s:
    Had[:,:,i] += (1-t)*Hi.full() 
    val, vec = np.linalg.eig(Had[:,:,i])

    
    sort_perm = val.argsort()

    val.sort()     # <-- This sorts the list in place.
    vec = vec[:, sort_perm]


    n[:,:,i] += vec # eigen_i = vec[:,i]
    E[i,:] += val
    i+=1
    
for i in range(len(s)-1):
    dt = s[i+1]-s[i]
    dx = n[:,:,i+1] - n[:,:,i]
    dn[:,:,i] += dx/dt

nn = np.zeros([2**l,2**l,len(s)],dtype=complex)
Enn = np.zeros([2**l,2**l,len(s)],dtype=complex)
dnn = np.zeros([2**l,2**l,len(s)],dtype=complex)
theta = np.zeros([len(s)],dtype=complex)
thtnn = np.zeros([2**l,2**l,len(s)],dtype=complex)
for i in range(len(s)-1):
    for j in range(2**l):
        Enn[:,:,i] += E[i,j]*np.outer(n[:,j,i],np.conjugate(n[:,j,i]))
        dnn[:,:,i] += np.outer(dn[:,j,i],np.conjugate(n[:,j,i]))
    
        theta[i] += np.dot(np.conjugate(n[:,j,i]),dn[:,j,i])
        nn[:,:,i] += np.outer(n[:,j,i],np.conjugate(n[:,j,i]))
        thtnn[:,:,i]+= (theta[i])*nn[:,:,i]


Hsad = Enn +1j*(dnn -thtnn)

In the analytical approach I've used SymPy. The problem here is to convert the symbolic expression I've founded to some that are possible to use as input to QuTip or odeintw.
but I have no idea on how to use it to solve my problem or if its possible. Can some one help me please?
 A: One option is to approximate your time-dependent Hamiltonian by a piecewise-constant Hamiltonian. Divide your time interval up into finite timesteps, and call odeintw or mesolve once for each step. The initial condition for step $i$ is $\rho$ at the end of step $i-1$.
A: You can solve your problem with qutip meSolve. Read this link to know how to handle time-dependents ones.
http://qutip.org/docs/latest/guide/dynamics/dynamics-time.html
Read part 2 & 3  :
1-String based: The Hamiltonian and/or collapse operators are expressed as a list of [qobj, string] pairs, where the time-dependent coefficients are represented as strings. The resulting Hamiltonian is then compiled into C code using Cython and executed.
2   -Array Based: The Hamiltonian and/or collapse operators are expressed as a list of [qobj, np.array] pairs. The arrays are 1 dimensional and dtype are complex or float. They must contain one value for each time in the tlist given to the solver. Cubic spline interpolation will be used between the given times.
