Why is this qutip code not working?

Question

I have been coding to calculate the following equation:

Where rho(t) is calculated by solving a quantum mechanics equation. <3| is a row vector with 1 at the 3rd element. |3> is a column vector with 1 at the 3rd element.

To solve the quantum mechanics equation (for the physics experts, I am solving the Lindblad Equation), I am using the qutip module.

The problem is, whenever I do the integration (using the scipy.integrate.quad function), I get a number like -0.842371561579 which is not what I expected (I am hoping for a number like 0.99). Also, I get a warning:

/usr/local/lib/python2.7/dist-packages/scipy/integrate/quadpack.py:352: IntegrationWarning: The integral is probably divergent, or slowly convergent.
  warnings.warn(msg, IntegrationWarning) 

I tried playing around with the values of the constants in the gamma constant. However, it didn't change any of the output. I have noticed it has to do with calculating rho(t), but that's all I could discover.

What is wrong with the code?

from qutip import *
import numpy as np
import scipy
from scipy.constants import *

hamiltonian = np.array([[215, -104.1, 5.1, -4.3  ,4.7,-15.1 ,-7.8 ],
[-104.1,  220.0, 32.6 ,7.1, 5.4, 8.3, 0.8],
      [ 5.1, 32.6, 0.0, -46.8, 1.0 , -8.1, 5.1],
     [-4.3, 7.1, -46.8, 125.0, -70.7, -14.7, -61.5],
       [ 4.7, 5.4, 1.0, -70.7, 450.0, 89.7, -2.5],
    [-15.1, 8.3, -8.1, -14.7, 89.7, 330.0, 32.7],
     [-7.8, 0.8, 5.1, -61.5,  -2.5, 32.7,  280.0]])
H=Qobj(hamiltonian) #This just converts hamiltonian into the qutip format

ground=H.groundstate()[1]

rho0 = ground*ground.dag()
gamma=(2*pi)*(296*0.695)*(35/150)
print gamma

#collapse operators:
L1 = Qobj(np.array([[1,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L2 = Qobj(np.array([[0,0,0,0,0,0,0],[0,1,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L3 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,1,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L4 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,1,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L5 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,1,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L6 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,1,0],[0,0,0,0,0,0,0]]))
L7 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,1]]))
#Since our gamma variable cannot be directly applied onto the Lindblad operator, we must multiply it with the collapse operators:

L1=gamma*L1
L2=gamma*L2
L3=gamma*L3
L4=gamma*L4
L5=gamma*L5
L6=gamma*L6
L7=gamma*L7

options = Options(nsteps=10000000, atol=1e-20)

#This code is for function-based integration
def integratefunc(x):
    results = mesolve(H, rho0, [x], [L1,L2,L3,L4,L5,L6,L7],[], options=options) #this function find rho(t) and evaluates it at t=x
    bra3= [[0,0,1,0,0,0,0]] #<3|
    bra3q=Qobj(bra3) #convert to qutip format
    ket3= [[0],[0],[1],[0],[0],[0],[0]] #|3>
    ket3q=Qobj(ket3) #convert to qutip format
    tempq=bra3q*(results.states[0]*ket3q) #multiply <3|rho(x)|3>
    y=tempq.data.data[0].real #gets the real part of the number, since the imaginary part is zero
    return y #returns the value

efficiency, error = scipy.integrate.quad(integratefunc, 0, np.inf) #does the integration to calculate the efficiency and the error
print efficiency, error #prints the final results