def _sesolve_list_td(H_func, psi0, tlist, e_ops, args, opt, progress_bar):
    """!
    Evolve the wave function using an ODE solver with time-dependent
    Hamiltonian.
    """

    if debug:
        print(inspect.stack()[0][3])

    if not isket(psi0):
        raise TypeError("psi0 must be a ket")

    #
    # configure time-dependent terms and setup ODE solver
    #
    if len(H_func) != 2:
        raise TypeError('Time-dependent Hamiltonian list must have two terms.')
    if (not isinstance(H_func[0], (list, np.ndarray))) or \
       (len(H_func[0]) <= 1):
        raise TypeError('Time-dependent Hamiltonians must be a list with two '
                        + 'or more terms')
    if (not isinstance(H_func[1], (list, np.ndarray))) or \
       (len(H_func[1]) != (len(H_func[0]) - 1)):
        raise TypeError('Time-dependent coefficients must be list with ' +
                        'length N-1 where N is the number of ' +
                        'Hamiltonian terms.')
    tflag = 1
    if opt.rhs_reuse and odeconfig.tdfunc is None:
        print("No previous time-dependent RHS found.")
        print("Generating one for you...")
        rhs_generate(H_func, args)
    lenh = len(H_func[0])
    if opt.tidy:
        H_func[0] = [(H_func[0][k]).tidyup() for k in range(lenh)]
    # create data arrays for time-dependent RHS function
    Hdata = [-1.0j * H_func[0][k].data.data for k in range(lenh)]
    Hinds = [H_func[0][k].data.indices for k in range(lenh)]
    Hptrs = [H_func[0][k].data.indptr for k in range(lenh)]
    # setup ode args string
    string = ""
    for k in range(lenh):
        string += ("Hdata[" + str(k) + "],Hinds[" + str(k) +
                   "],Hptrs[" + str(k) + "],")

    if args:
        td_consts = args.items()
        for elem in td_consts:
            string += str(elem[1])
            if elem != td_consts[-1]:
                string += (",")

    # run code generator
    if not opt.rhs_reuse or odeconfig.tdfunc is None:
        if opt.rhs_filename is None:
            odeconfig.tdname = "rhs" + str(odeconfig.cgen_num)
        else:
            odeconfig.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args,
                       odeconfig=odeconfig)
        cgen.generate(odeconfig.tdname + ".pyx")

        code = compile('from ' + odeconfig.tdname + ' import cy_td_ode_rhs',
                       '<string>', 'exec')
        exec(code, globals())
        odeconfig.tdfunc = cy_td_ode_rhs
    #
    # setup integrator
    #
    initial_vector = psi0.full().ravel()
    r = scipy.integrate.ode(odeconfig.tdfunc)
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
                     first_step=opt.first_step, min_step=opt.min_step,
                     max_step=opt.max_step)
    r.set_initial_value(initial_vector, tlist[0])
    code = compile('r.set_f_params(' + string + ')', '<string>', 'exec')
    exec(code)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, psi0, tlist, e_ops, opt, progress_bar,
                              norm, dims=psi0.dims)

def _mesolve_list_td(H_func, rho0, tlist, c_op_list, e_ops, args, opt, progress_bar):
    """
    Evolve the density matrix using an ODE solver with time dependent
    Hamiltonian.
    """

    if debug:
        print(inspect.stack()[0][3])

    #
    # check initial state
    #
    if isket(rho0):
        # if initial state is a ket and no collapse operator where given,
        # fall back on the unitary schrodinger equation solver
        if len(c_op_list) == 0:
            return _sesolve_list_td(H_func, rho0, tlist, e_ops, args, opt, progress_bar)

        # Got a wave function as initial state: convert to density matrix.
        rho0 = ket2dm(rho0)

    #
    # construct liouvillian
    #
    if len(H_func) != 2:
        raise TypeError("Time-dependent Hamiltonian list must have two terms.")
    if not isinstance(H_func[0], (list, np.ndarray)) or len(H_func[0]) <= 1:
        raise TypeError("Time-dependent Hamiltonians must be a list " + "with two or more terms")
    if (not isinstance(H_func[1], (list, np.ndarray))) or (len(H_func[1]) != (len(H_func[0]) - 1)):
        raise TypeError(
            "Time-dependent coefficients must be list with "
            + "length N-1 where N is the number of "
            + "Hamiltonian terms."
        )

    if opt.rhs_reuse and config.tdfunc is None:
        rhs_generate(H_func, args)

    lenh = len(H_func[0])
    if opt.tidy:
        H_func[0] = [(H_func[0][k]).tidyup() for k in range(lenh)]
    L_func = [[liouvillian(H_func[0][0], c_op_list)], H_func[1]]
    for m in range(1, lenh):
        L_func[0].append(liouvillian(H_func[0][m], []))

    # create data arrays for time-dependent RHS function
    Ldata = [L_func[0][k].data.data for k in range(lenh)]
    Linds = [L_func[0][k].data.indices for k in range(lenh)]
    Lptrs = [L_func[0][k].data.indptr for k in range(lenh)]
    # setup ode args string
    string = ""
    for k in range(lenh):
        string += "Ldata[%d], Linds[%d], Lptrs[%d]," % (k, k, k)

    if args:
        td_consts = args.items()
        for elem in td_consts:
            string += str(elem[1])
            if elem != td_consts[-1]:
                string += ","

    # run code generator
    if not opt.rhs_reuse or config.tdfunc is None:
        if opt.rhs_filename is None:
            config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)
        else:
            config.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args, config=config)
        cgen.generate(config.tdname + ".pyx")

        code = compile("from " + config.tdname + " import cy_td_ode_rhs", "<string>", "exec")
        exec(code, globals())
        config.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel()
    r = scipy.integrate.ode(config.tdfunc)
    r.set_integrator(
        "zvode",
        method=opt.method,
        order=opt.order,
        atol=opt.atol,
        rtol=opt.rtol,
        nsteps=opt.nsteps,
        first_step=opt.first_step,
        min_step=opt.min_step,
        max_step=opt.max_step,
    )
    r.set_initial_value(initial_vector, tlist[0])
    code = compile("r.set_f_params(" + string + ")", "<string>", "exec")
    exec(code)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)

def _sesolve_list_str_td(H_list, psi0, tlist, e_ops, args, opt,
                         progress_bar):
    """
    Internal function for solving the master equation. See mesolve for usage.
    """

    if debug:
        print(inspect.stack()[0][3])

    #
    # check initial state: must be a density matrix
    #
    if not isket(psi0):
        raise TypeError("The unitary solver requires a ket as initial state")

    #
    # construct liouvillian
    #
    Ldata = []
    Linds = []
    Lptrs = []
    Lcoeff = []

    # loop over all hamiltonian terms, convert to superoperator form and
    # add the data of sparse matrix representation to h_coeff
    for h_spec in H_list:

        if isinstance(h_spec, Qobj):
            h = h_spec
            h_coeff = "1.0"

        elif isinstance(h_spec, list):
            h = h_spec[0]
            h_coeff = h_spec[1]

        else:
            raise TypeError("Incorrect specification of time-dependent " +
                            "Hamiltonian (expected string format)")

        L = -1j * h

        Ldata.append(L.data.data)
        Linds.append(L.data.indices)
        Lptrs.append(L.data.indptr)
        Lcoeff.append(h_coeff)

    # the total number of liouvillian terms (hamiltonian terms +
    # collapse operators)
    n_L_terms = len(Ldata)

    #
    # setup ode args string: we expand the list Ldata, Linds and Lptrs into
    # and explicit list of parameters
    #
    string_list = []
    for k in range(n_L_terms):
        string_list.append("Ldata[%d], Linds[%d], Lptrs[%d]" % (k, k, k))
    for name, value in args.items():
        string_list.append(str(value))
    parameter_string = ",".join(string_list)

    #
    # generate and compile new cython code if necessary
    #
    if not opt.rhs_reuse or odeconfig.tdfunc is None:
        if opt.rhs_filename is None:
            odeconfig.tdname = "rhs" + str(odeconfig.cgen_num)
        else:
            odeconfig.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args,
                       odeconfig=odeconfig)
        cgen.generate(odeconfig.tdname + ".pyx")

        code = compile('from ' + odeconfig.tdname + ' import cy_td_ode_rhs',
                       '<string>', 'exec')
        exec(code, globals())
        odeconfig.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
    initial_vector = psi0.full().ravel()
    r = scipy.integrate.ode(odeconfig.tdfunc)
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
                     first_step=opt.first_step, min_step=opt.min_step,
                     max_step=opt.max_step)
    r.set_initial_value(initial_vector, tlist[0])
    code = compile('r.set_f_params(' + parameter_string + ')',
                   '<string>', 'exec')
    exec(code)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, psi0, tlist, e_ops, opt, progress_bar,
                              norm, dims=psi0.dims)

#.........这里部分代码省略.........
        col_spmv_code = ("state = _cy_col_spmv_func(j, ODE.t, " +
                         "config.c_ops_data[j], config.c_ops_ind[j], " +
                         "config.c_ops_ptr[j], ODE.y")
        col_expect_code = ("for i in config.c_td_inds: " +
                           "n_dp.append(_cy_col_expect_func(i, ODE.t, " +
                           "config.n_ops_data[i], " +
                           "config.n_ops_ind[i], " +
                           "config.n_ops_ptr[i], ODE.y")
        for kk in range(len(config.c_args)):
            col_spmv_code += ",config.c_args[" + str(kk) + "]"
            col_expect_code += ",config.c_args[" + str(kk) + "]"
        col_spmv_code += ")"
        col_expect_code += "))"

        config.col_spmv_code = col_spmv_code
        config.col_expect_code = col_expect_code

        # setup ode args string
        config.string = ""
        data_range = range(len(config.h_data))
        for k in data_range:
            config.string += ("config.h_data[" + str(k) +
                              "], config.h_ind[" + str(k) +
                              "], config.h_ptr[" + str(k) + "]")
            if k != data_range[-1]:
                config.string += ","
        # attach args to ode args string
        if len(config.c_args) > 0:
            for kk in range(len(config.c_args)):
                config.string += "," + "config.c_args[" + str(kk) + "]"

        name = "rhs" + str(os.getpid()) + str(config.cgen_num)
        config.tdname = name
        cgen = Codegen(H_inds, H_tdterms, config.h_td_inds, args,
                       C_inds, C_tdterms, config.c_td_inds, type='mc',
                       config=config)
        cgen.generate(name + ".pyx")

    elif config.tflag in [2, 20, 22]:
        # PYTHON LIST-FUNCTION BASED TIME-DEPENDENCE
        # ------------------------------------------

        # take care of Hamiltonian
        if config.tflag == 2:
            # constant Hamiltonian, at least one function based collapse
            # operators
            H_inds = np.array([0])
            H_tdterms = 0
            len_h = 1
        else:
            # function based Hamiltonian
            H_inds = np.arange(len(H))
            H_td_inds = np.array(h_stuff[1])
            H_const_inds = np.setdiff1d(H_inds, H_td_inds)
            config.h_funcs = np.array([H[k][1] for k in H_td_inds])
            config.h_func_args = args
            Htd = np.array([H[k][0] for k in H_td_inds], dtype=object)
            config.h_td_inds = np.arange(len(Htd))
            H = np.sum(H[k] for k in H_const_inds)

        # take care of collapse operators
        C_inds = np.arange(config.c_num)
        # find inds of time-dependent terms
        C_td_inds = np.array(c_stuff[1])
        # find inds of constant terms
        C_const_inds = np.setdiff1d(C_inds, C_td_inds)

#.........这里部分代码省略.........
                Lconst += c
            else:
                raise TypeError(
                    "Incorrect specification of time-dependent "
                    + "Liouvillian (expected operator or "
                    + "superoperator)"
                )

        elif isinstance(c_spec, list):
            c = c_spec[0]
            c_coeff = c_spec[1]

            if isoper(c):
                cdc = c.dag() * c
                L = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
                c_coeff = "(" + c_coeff + ")**2"
            elif issuper(c):
                L = c
            else:
                raise TypeError(
                    "Incorrect specification of time-dependent "
                    + "Liouvillian (expected operator or "
                    + "superoperator)"
                )

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append(c_coeff)

        else:
            raise TypeError(
                "Incorrect specification of time-dependent " + "collapse operators (expected string format)"
            )

    # add the constant part of the lagrangian
    if Lconst != 0:
        Ldata.append(Lconst.data.data)
        Linds.append(Lconst.data.indices)
        Lptrs.append(Lconst.data.indptr)
        Lcoeff.append("1.0")

    # the total number of liouvillian terms (hamiltonian terms +
    # collapse operators)
    n_L_terms = len(Ldata)

    #
    # setup ode args string: we expand the list Ldata, Linds and Lptrs into
    # and explicit list of parameters
    #
    string_list = []
    for k in range(n_L_terms):
        string_list.append("Ldata[%d], Linds[%d], Lptrs[%d]" % (k, k, k))
    for name, value in args.items():
        if isinstance(value, np.ndarray):
            string_list.append(name)
        else:
            string_list.append(str(value))
    parameter_string = ",".join(string_list)

    #
    # generate and compile new cython code if necessary
    #
    if not opt.rhs_reuse or config.tdfunc is None:
        if opt.rhs_filename is None:
            config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)
        else:
            config.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args, config=config)
        cgen.generate(config.tdname + ".pyx")

        code = compile("from " + config.tdname + " import cy_td_ode_rhs", "<string>", "exec")
        exec(code, globals())
        config.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel()
    r = scipy.integrate.ode(config.tdfunc)
    r.set_integrator(
        "zvode",
        method=opt.method,
        order=opt.order,
        atol=opt.atol,
        rtol=opt.rtol,
        nsteps=opt.nsteps,
        first_step=opt.first_step,
        min_step=opt.min_step,
        max_step=opt.max_step,
    )
    r.set_initial_value(initial_vector, tlist[0])
    code = compile("r.set_f_params(" + parameter_string + ")", "<string>", "exec")

    exec(code, locals(), args)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)

#.........这里部分代码省略.........
            if not isinstance(h, Qobj):
                raise TypeError(msg + "expected Qobj")

            if h.isoper:
                Lconst += -1j * (spre(h) - spost(h))
            elif h.issuper:
                Lconst += h
            else:
                raise TypeError(msg + "expected operator or superoperator")

        elif isinstance(h_spec, list):
            h = h_spec[0]
            h_coeff = h_spec[1]

            if not isinstance(h, Qobj):
                raise TypeError(msg + "expected Qobj")

            if h.isoper:
                L = -1j * (spre(h) - spost(h))
            elif h.issuper:
                L = h
            else:
                raise TypeError(msg + "expected operator or superoperator")

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append(h_coeff)

        else:
            raise TypeError(msg + "expected string format")

    # loop over all collapse operators
    for c_spec in c_ops:
        if isinstance(c_spec, Qobj):
            c = c_spec

            if not isinstance(c, Qobj):
                raise TypeError(msg + "expected Qobj")

            if c.isoper:
                cdc = c.dag() * c
                Lconst += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) \
                                                   - 0.5 * spost(cdc)
            elif c.issuper:
                Lconst += c
            else:
                raise TypeError(msg + "expected operator or superoperator")

        elif isinstance(c_spec, list):
            c = c_spec[0]
            c_coeff = c_spec[1]

            if not isinstance(c, Qobj):
                raise TypeError(msg + "expected Qobj")

            if c.isoper:
                cdc = c.dag() * c
                L = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) \
                                             - 0.5 * spost(cdc)
                c_coeff = "(" + c_coeff + ")**2"
            elif c.issuper:
                L = c
            else:
                raise TypeError(msg + "expected operator or superoperator")

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append(c_coeff)

        else:
            raise TypeError(msg + "expected string format")

    # add the constant part of the lagrangian
    if Lconst != 0:
        Ldata.append(Lconst.data.data)
        Linds.append(Lconst.data.indices)
        Lptrs.append(Lconst.data.indptr)
        Lcoeff.append("1.0")

    # the total number of liouvillian terms (hamiltonian terms + collapse
    # operators)
    n_L_terms = len(Ldata)

    cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args,
                   config=config)
    cgen.generate(config.tdname + ".pyx")

    code = compile('from ' + config.tdname +
                   ' import cy_td_ode_rhs', '<string>', 'exec')
    exec(code, globals())

    config.tdfunc = cy_td_ode_rhs

    if cleanup:
        try:
            os.remove(config.tdname + ".pyx")
        except:
            pass

def _sesolve_list_str_td(H_list, psi0, tlist, e_ops, args, opt,
                         progress_bar):
    """
    Internal function for solving the master equation. See mesolve for usage.
    """

    if debug:
        print(inspect.stack()[0][3])

    if psi0.isket:
        oper_evo = False
    elif psi0.isunitary:
        oper_evo = True
    else:
        raise TypeError("The unitary solver requires psi0 to be"
                        " a ket as initial state"
                        " or a unitary as initial operator.")

    #
    # construct dynamics generator
    #
    Ldata = []
    Linds = []
    Lptrs = []
    Lcoeff = []
    Lobj = []

    # loop over all hamiltonian terms, convert to superoperator form and
    # add the data of sparse matrix representation to h_coeff
    for h_spec in H_list:

        if isinstance(h_spec, Qobj):
            h = h_spec
            h_coeff = "1.0"

        elif isinstance(h_spec, list):
            h = h_spec[0]
            h_coeff = h_spec[1]

        else:
            raise TypeError("Incorrect specification of time-dependent " +
                            "Hamiltonian (expected string format)")

        L = -1j * h
        Ldata.append(L.data.data)
        Linds.append(L.data.indices)
        Lptrs.append(L.data.indptr)
        if isinstance(h_coeff, Cubic_Spline):
            Lobj.append(h_coeff.coeffs)
        Lcoeff.append(h_coeff)

    # the total number of Hamiltonian terms
    n_L_terms = len(Ldata)

    # Check which components should use OPENMP
    omp_components = None
    if qset.has_openmp:
        if opt.use_openmp:
            omp_components = openmp_components(Lptrs)

    #
    # setup ode args string: we expand the list Ldata, Linds and Lptrs into
    # and explicit list of parameters
    #
    string_list = []
    for k in range(n_L_terms):
        string_list.append("Ldata[%d], Linds[%d], Lptrs[%d]" % (k, k, k))
    # Add object terms to end of ode args string
    for k in range(len(Lobj)):
        string_list.append("Lobj[%d]" % k)

    for name, value in args.items():
        if isinstance(value, np.ndarray):
            string_list.append(name)
        else:
            string_list.append(str(value))
    parameter_string = ",".join(string_list)

    #
    # generate and compile new cython code if necessary
    #
    if not opt.rhs_reuse or config.tdfunc is None:
        if opt.rhs_filename is None:
            config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)
        else:
            config.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args,
                       config=config, use_openmp=opt.use_openmp,
                       omp_components=omp_components,
                       omp_threads=opt.openmp_threads)
        cgen.generate(config.tdname + ".pyx")

        code = compile('from ' + config.tdname + ' import cy_td_ode_rhs',
                       '<string>', 'exec')
        exec(code, globals())
        config.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
#.........这里部分代码省略.........

#.........这里部分代码省略.........
                                "superoperator)")

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            #Lcoeff.append(c_coeff)

        else:
            raise TypeError("Incorrect specification of time-dependent " +
                            "collapse operators (expected string format)")
    
    
    #prepend the constant part of the liouvillian
    if Lconst != 0:
       Ldata = [Lconst.data.data]+Ldata
       Linds = [Lconst.data.indices]+Linds
       Lptrs = [Lconst.data.indptr]+Lptrs
       Lcoeff = ["1.0"]+Lcoeff
       
    else:
        me_cops_obj_flags = [kk-1 for kk in me_cops_obj_flags]
    # the total number of liouvillian terms (hamiltonian terms +
    # collapse operators)
    n_L_terms = len(Ldata)
    n_td_cops = len(me_cops_obj)
    
    # Check which components should use OPENMP
    omp_components = None
    if qset.has_openmp:
        if opt.use_openmp:
            omp_components = openmp_components(Lptrs)

    #
    # setup ode args string: we expand the list Ldata, Linds and Lptrs into
    # and explicit list of parameters
    #
    string_list = []
    for k in range(n_L_terms):
        string_list.append("Ldata[%d], Linds[%d], Lptrs[%d]" % (k, k, k))
    
    # Add H object terms to ode args string
    for k in range(len(Lobj)):
        string_list.append("Lobj[%d]" % k)
        
    # Add cop object terms to end of ode args string
    for k in range(len(me_cops_obj)):
        string_list.append("me_cops_obj[%d]" % k)    
    
    for name, value in args.items():
        if isinstance(value, np.ndarray):
            string_list.append(name)
        else:
            string_list.append(str(value))
    parameter_string = ",".join(string_list)
    
    #
    # generate and compile new cython code if necessary
    #
    if not opt.rhs_reuse or config.tdfunc is None:
        if opt.rhs_filename is None:
            config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)
        else:
            config.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=len(Lcoeff), h_tdterms=Lcoeff, 
                       c_td_splines=me_cops_coeff, 
                       c_td_spline_flags=me_cops_obj_flags, args=args,
                       config=config, use_openmp=opt.use_openmp,
                       omp_components=omp_components,
                       omp_threads=opt.openmp_threads)
        cgen.generate(config.tdname + ".pyx")

        code = compile('from ' + config.tdname + ' import cy_td_ode_rhs',
                       '<string>', 'exec')
        exec(code, globals())
        config.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel('F')
    if issuper(rho0):
        r = scipy.integrate.ode(_td_ode_rhs_super)
        code = compile('r.set_f_params([' + parameter_string + '])',
                       '<string>', 'exec')
    else:
        r = scipy.integrate.ode(config.tdfunc)
        code = compile('r.set_f_params(' + parameter_string + ')',
                       '<string>', 'exec')
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
                     first_step=opt.first_step, min_step=opt.min_step,
                     max_step=opt.max_step)
    r.set_initial_value(initial_vector, tlist[0])

    exec(code, locals(), args)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)

def rhs_generate(H, c_ops, args={}, options=Odeoptions(), name=None):
    """
    Generates the Cython functions needed for solving the dynamics of a
    given system using the mesolve function inside a parfor loop.

    Parameters
    ----------
    H : qobj
        System Hamiltonian.
    c_ops : list
        ``list`` of collapse operators.
    args : dict
        Arguments for time-dependent Hamiltonian and collapse operator terms.
    options : Odeoptions
        Instance of ODE solver options.
    name: str
        Name of generated RHS

    Notes
    -----
    Using this function with any solver other than the mesolve function
    will result in an error.

    """
    odeconfig.reset()
    odeconfig.options = options

    if name:
        odeconfig.tdname = name
    else:
        odeconfig.tdname = "rhs" + str(odeconfig.cgen_num)

    Lconst = 0

    Ldata = []
    Linds = []
    Lptrs = []
    Lcoeff = []

    # loop over all hamiltonian terms, convert to superoperator form and
    # add the data of sparse matrix represenation to
    for h_spec in H:
        if isinstance(h_spec, Qobj):
            h = h_spec
            Lconst += -1j * (spre(h) - spost(h))

        elif isinstance(h_spec, list):
            h = h_spec[0]
            h_coeff = h_spec[1]

            L = -1j * (spre(h) - spost(h))

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append(h_coeff)

        else:
            raise TypeError("Incorrect specification of time-dependent " + "Hamiltonian (expected string format)")

    # loop over all collapse operators
    for c_spec in c_ops:
        if isinstance(c_spec, Qobj):
            c = c_spec
            cdc = c.dag() * c
            Lconst += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)

        elif isinstance(c_spec, list):
            c = c_spec[0]
            c_coeff = c_spec[1]

            cdc = c.dag() * c
            L = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append("(" + c_coeff + ")**2")

        else:
            raise TypeError(
                "Incorrect specification of time-dependent " + "collapse operators (expected string format)"
            )

    # add the constant part of the lagrangian
    if Lconst != 0:
        Ldata.append(Lconst.data.data)
        Linds.append(Lconst.data.indices)
        Lptrs.append(Lconst.data.indptr)
        Lcoeff.append("1.0")

    # the total number of liouvillian terms (hamiltonian terms + collapse
    # operators)
    n_L_terms = len(Ldata)

    cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args, odeconfig=odeconfig)
    cgen.generate(odeconfig.tdname + ".pyx")

    code = compile("from " + odeconfig.tdname + " import cy_td_ode_rhs", "<string>", "exec")
    exec(code, globals())
#.........这里部分代码省略.........

def _mesolve_list_td(H_func, rho0, tlist, c_op_list, e_ops, args, opt,
                     progress_bar):
    """!
    Evolve the density matrix using an ODE solver with time dependent
    Hamiltonian.
    """

    if debug:
        print(inspect.stack()[0][3])

    #
    # check initial state
    #
    if isket(rho0):
        # if initial state is a ket and no collapse operator where given,
        # fall back on the unitary schrodinger equation solver
        if len(c_op_list) == 0:
            return _sesolve_list_td(H_func, rho0, tlist, e_ops, args, opt)

        # Got a wave function as initial state: convert to density matrix.
        rho0 = ket2dm(rho0)

    #
    # construct liouvillian
    #
    if len(H_func) != 2:
        raise TypeError('Time-dependent Hamiltonian list must have two terms.')
    if not isinstance(H_func[0], (list, np.ndarray)) or len(H_func[0]) <= 1:
        raise TypeError('Time-dependent Hamiltonians must be a list ' +
                        'with two or more terms')
    if (not isinstance(H_func[1], (list, np.ndarray))) or \
       (len(H_func[1]) != (len(H_func[0]) - 1)):
        raise TypeError('Time-dependent coefficients must be list with ' +
                        'length N-1 where N is the number of ' +
                        'Hamiltonian terms.')

    if opt.rhs_reuse and odeconfig.tdfunc is None:
        rhs_generate(H_func, args)

    lenh = len(H_func[0])
    if opt.tidy:
        H_func[0] = [(H_func[0][k]).tidyup() for k in range(lenh)]
    L_func = [[liouvillian_fast(H_func[0][0], c_op_list)], H_func[1]]
    for m in range(1, lenh):
        L_func[0].append(liouvillian_fast(H_func[0][m], []))

    # create data arrays for time-dependent RHS function
    Ldata = [L_func[0][k].data.data for k in range(lenh)]
    Linds = [L_func[0][k].data.indices for k in range(lenh)]
    Lptrs = [L_func[0][k].data.indptr for k in range(lenh)]
    # setup ode args string
    string = ""
    for k in range(lenh):
        string += ("Ldata[%d], Linds[%d], Lptrs[%d]," % (k, k, k))

    if args:
        for name, value in args.items():
            if isinstance(value, np.ndarray):
                globals()['var_%s'%name] = value
                string += 'var_%s,'%name
            else:
                string += str(value) + ','

    # run code generator
    if not opt.rhs_reuse or odeconfig.tdfunc is None:
        if opt.rhs_filename is None:
            odeconfig.tdname = "rhs" + str(odeconfig.cgen_num)
        else:
            odeconfig.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args,
                       odeconfig=odeconfig)
        cgen.generate(odeconfig.tdname + ".pyx")

        code = compile('from ' + odeconfig.tdname + ' import cy_td_ode_rhs',
                       '<string>', 'exec')
        exec(code, globals())
        odeconfig.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel()
    r = scipy.integrate.ode(odeconfig.tdfunc)
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
                     first_step=opt.first_step, min_step=opt.min_step,
                     max_step=opt.max_step)
    r.set_initial_value(initial_vector, tlist[0])
    code = compile('r.set_f_params(' + string + ')', '<string>', 'exec')
    exec(code)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)