Input output

import numpy as np
import matplotlib.pyplot as plt
import qutip as qt
from scipy.integrate import quad

def u(t):
    sigma = 2/gamma
    return np.exp(-t**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))

def v(t):
    sigma = 2/gamma
    return np.exp(-t**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))

def integrandu(t):
    return np.abs(u(t))**2

def integrandv(t):
    return np.abs(v(t))**2


def gu(ti):
    return np.conj(u(ti))/np.sqrt(1-quad(integrandu, t[0], ti)[0])

def gv(ti):
    return np.conj(v(ti))/np.sqrt(quad(integrandv, t[0], ti)[0])

gamma = 2*np.pi*1
wc = 2*np.pi*484115/242

t = np.linspace(0, 12*1/gamma, 10000)


psi0 = qt.tensor(qt.basis(2, 1), qt.basis(2, 0), qt.basis(2, 0))
rho0 = psi0*psi0.dag()

au = qt.tensor(qt.destroy(2), qt.qeye(2), qt.qeye(2))
c = qt.tensor(qt.qeye(2), qt.destroy(2), qt.qeye(2))
av = qt.tensor(qt.qeye(2), qt.qeye(2), qt.destroy(2))
L0 = list(np.zeros(np.shape(t)))
for k, ti in enumerate(t):
    L0[k] = np.sqrt(gamma)*c+np.conj(gu(ti))*au+np.conj(gv(ti))*av


Hs = wc*c.dag()*c
H = list(np.zeros(np.shape(t)))
for k, ti in enumerate(t):
    cascade = (np.sqrt(gamma)*gu(ti)*au.dag()*c+
                np.sqrt(gamma)*gv(ti)*c.dag()*av+
                gv(ti)*np.conj(gv(ti))*au.dag()*av)
    cascade = (1j/2)*(cascade-cascade.dag())
    H[k] = Hs+cascade


c_ops = [L0]
e_ops = [au.dag()*au, c.dag()*c, av.dag()*av]

output = qt.mesolve(H, rho0, t, L0, e_ops)

plt.plot(t, output.expect[0])
plt.show()