<p>From a programming point of view, the problem appears to be your value of <code>gamma</code> and therefore the size of your collapse operators.  Print out <code>gamma</code> - it is of the order <code>10**25</code> - this seems to be what is preventing the solver from converging.</p>

<p>Just for testing (I'm an engineer, not a quantum physicist...), I put in a smaller value of <code>gamma</code> (e.g. 0.1), the solver seems to work and gives apparently reasonable output in <code>results.states</code></p>

<p>I don't quite understand your <code>gamma</code> - it <em>seems</em> to have units of cm<sup>-1</sup>s<sup>-2</sup> as you have set it up.  I wonder if you only want to divide by <code>hbar</code> once, maybe.  As I say, I'm not a quantum physicist, so I'm only guessing here based on what makes the programming hang together combined with a bit of dimensional analysis.</p>

<p><strong><em>EDIT</em></strong></p>

<p>OP indicates in comments that the wrong order of magnitude / units for <code>gamma</code> does seem to be the programming issue (i.e. preventing numerical calculus from converging), but isn't totally clear on how to calculate gamma.  At this stage, it may be worth posting a question at either <a href="http://physics.stackexchange.com">http://physics.stackexchange.com</a> or <a href="http://math.stackexchange.com">http://math.stackexchange.com</a> about that - referencing this one for context if necessary.</p>

<p><strong><em>EDIT 2</em></strong></p>

<p>I note you asked this <a href="https://physics.stackexchange.com/questions/187734/messed-up-units">related question on the physics site</a>.  This makes it clear where <a href="http://arxiv.org/pdf/0807.0929v2.pdf" rel="nofollow noreferrer">the expression for gamma comes from</a> and thereby clarifies that the constant terms presented as simply <code>30</code> and <code>150</code> in this question actually have units (Energy and frequency respectively).  This changes the dimensional analysis - the units of gamma are s<sup>-1</sup> or, with appropriate conversion, cm<sup>-1</sup>.</p>

<p>It also shows the value you mention in comments - 300 cm<sup>-1</sup>.</p>
Integration not successful in Python QuTiP
<p>I have been trying to use QuTiP to solve a quantum mechanics matrix differential equation (a Lindblad equation). Here is the code:</p>

<pre><code>from qutip import *
from matplotlib import *
import numpy as np

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)
ground=Qobj(np.array([[ 0.0863685 ],
  [ 0.17141713],
  [-0.91780802],
  [-0.33999268],
  [-0.04835763],
  [-0.01859027],
  [-0.05006013]]))

rho0 = ground*ground.dag()
from scipy.constants import *
ktuple=physical_constants['Boltzmann constant in eV/K']
k = ktuple[0]* 8065.6
htuple = physical_constants['Planck constant in eV s'] 
hbar = (htuple[0]* 8065.6)/(2*pi)
gamma=(2*pi)*((k*300)/hbar)*(35/(150*hbar))

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=100000)

results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
print results
</code></pre>

<p>This code is supposed to solve the following equation:</p>

<p><img src="https://i.stack.imgur.com/x2X0Z.gif" alt="Lindblad equation"></p>

<p>where L_i are matrices (in the list: [L1,L2,L3,L4,L5,L6,L7]), H is the hamiltonian, another matrix, <img src="https://i.stack.imgur.com/TeKxt.gif" alt="$\rho$"> is the density matrix, and <img src="https://i.stack.imgur.com/C2zjo.gif" alt="$\gamma$"> is a constant equal to <img src="https://i.stack.imgur.com/JaFwh.gif" alt="$2\pi  kT/\hbar*E_{R}/(\hbar\omega_{c})$"> where T is the temperature, k is the Boltzmann constant, and <img src="https://i.stack.imgur.com/BYYcG.gif" alt="$\hbar$ = $h/2\pi$">, where h is Planck's constant.</p>

<p>Every time I run the code, it gives me the following error:</p>

<pre><code>ZVODE--  At T (=R1) and step size H (=R2), the
       corrector convergence failed repeatedly
       or with abs(H) = HMIN
      In above,  R1 =  0.0000000000000D+00   R2 =  0.1202322246215D-36
/usr/local/lib/python2.7/dist-packages/scipy/integrate/_ode.py:853: UserWarning: zvode: Repeated convergence failures. (Perhaps bad Jacobian supplied or wrong choice of MF o
r tolerances.)
  'Unexpected istate=%s' % istate))
Traceback (most recent call last):
  File "lindbladqutip.py", line 48, in &lt;module&gt;
    results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 264, in mesolve
    progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 692, in _mesolve_const
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 866, in _generic_ode_solve
    raise Exception("ODE integration error: Try to increase "
Exception: ODE integration error: Try to increase the allowed number of substeps by increasing the nsteps parameter in the Options class.
</code></pre>

<p>After doing some debugging analysis, it seems like the first or second integration fails. The error tells me to increase the nsteps parameter, which I have tried. Even then it fails. Changing the list of times (the np.linspace function makes the list of times) also has no effect. </p>

<p>I am desperate to know what I can do to fix this error. Please comment if you all need more details.</p>

<p>Thanks for all your help!</p>