Eigenvalues calculated using Qutip breaks for specific values of coupling constant

Question

I am trying to create a plot of Eigenenergies vs coupling constant g for my Jaynes-Cummings (JC) hamiltonian. I have noticed that when I try this, the eigenenergies seem to "break" at points where I would expect there to be crossings between eigenenergies. I am new to Qutip, so this could very well be an error in my code, however I have not been able to spot it. Is there a way to fix this issue?

Here is a portion of my code:

def JC(N, g, omega, delta):

    sz = tensor(qeye(N), sigmaz())
    a = tensor(destroy(N), qeye(2))

    sp = tensor(qeye(N), sigmap())
    sm = tensor(qeye(N), sigmam())

    H = omega*a.dag()*a + 0.5*(delta + omega)*sz + g*(a*sp + a.dag()*sm)

    EigVal, EigVec = H.eigenstates()

    return EigVal

When I plot the i'th eigenenergies vs g (for omega = 1 and detuning = 0), I get

I have tried different ways of writing the operators, and tried using a different technique to calculate the eigenvalues of the hamiltonian but other problems arose. I have calculated this analytically and I expect a similar plot however without the zigzag features.