title: Excess work done on call (perhaps wrong Dfun type) with scipy.integrate.odeint
body: <p>I'm writing a code in Python that predicts the energy levels of Hydrogen which I will use as a template for research into quarkonium energy levels. I'm using the <code>scipy.integrate.odeint()</code> function to solve the Shroedinger equation and it works fine for the lower energy levels up to <code>n=6</code>. I don't expect I'll have much need to go beyond that, but odeint returns <code>Excess work done on this call (perhaps wrong Dfun type).</code> which only encourages me to extend what I can predict.</p>
  
  <p>The Shroedinger equation substitution I'm using is:</p>
  
  <pre><code>u'' - (l*(l+1)/r**2 - 2mu_e(E-V_emag(r))) * u = 0
  =&gt;
  u' = v
  v' = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*u
  </code></pre>
  
  <p>I'm then using <code>scipy.integrate.odeint()</code> on it and iterating through energies and using other functions I've defined that assess turning points and nodes in the result. The way I find the energy levels is finding the lowest possible E value where the number of turning points and nodes matches what it should; then incrementing <code>L</code> by 1 and finding the new ground energy, e.g. if <code>L=0</code> I'll find <code>n=1</code> energy and if <code>L=3</code>, I'll find the <code>n=2</code> energy.</p>
  
  <p>Once the code increments to <code>L=7</code> it doesn't return anything useful. The range of <code>r</code> has been extended but I've tried with keeping it the same to reduce the number of steps but to no avail. The code is self-taught so in my research I've read about Jacobians. I'm afraid I haven't worked out what they are or if I need one yet. Any ideas?</p>
  
  <pre><code>def v_emag(r):
      v = -alpha/r
      return v
  
  def s_e(y,r,l,E): #Shroedinger equation for electromagntism
      x = numpy.zeros_like(y)
      x[0] = y[1]
      x[1] = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*y[0]
      return x
  
  def i_s_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
      if stop is None:
          stop = ((l+1)*30-10)*bohr
      r = numpy.arange(start,stop,step)
      y = odeint(s_e,y0,r,args=(l,E))
      return y
  
  def inormalise_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
      if stop is None:
          stop = ((l+1)*30-10)*bohr
      r = numpy.arange(start,stop,step)
      f = i_s_e(l,E,start,stop,step)[:,0]
      f2 = f**2
      area = numpy.trapz(f2,x=r)
      return f/(numpy.sqrt(area))
  
  def inodes_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
      if stop is None:
          stop = ((l+1)*30-10)*bohr
      x = i_s_e(l,E,start,stop,step)
      r = numpy.arange(start,stop,step)
      k=0
      for i in range(len(r)-1):
          if x[i,0]*x[i+1,0] &lt; 0: #If u value times next u value &lt;0,
              k+=1               #crossing of u=0 has occured, therefore count node
      return k
  
  def iturns_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
      if stop is None:
          stop = ((l+1)*30-10)*bohr
      x = i_s_e(l,E,start,stop,step)
      r = numpy.arange(start,stop,step)
      k = 0
      for i in range(len(r)-1):
          if x[i,1]*x[i+1,1] &lt; 0: #If du/dr value times next du/dr value &lt;0,
              k=k+1               #crossing of du/dr=0, therefore a maximum/minimum
      return k
  
  l = 0
  while l &lt; 10:    #The ground state for a system with a non-zero angular momentum will
      E1 = -1.5e-08    #be the energy level of principle quantum number l-1, therefore
      E3 = 0           #by changing l, we can find n by searching for the ground state
      E2 = 0.5*(E1+E3)
      i = 0
      while i &lt; 40:
          N1 = inodes_e(l,E1)
          N2 = inodes_e(l,E2)
          N3 = inodes_e(l,E3)
          T1 = iturns_e(l,E1)
          T2 = iturns_e(l,E2)
          T3 = iturns_e(l,E3)
          if N1 != N2:# and T1 != T2: #Looks in lower half first, therefore will tend to ground state
              E3 = E2
              E2 = 0.5*(E1+E3)
          elif N2 != N3:# and T2 != T3:
              E1 = E2
              E2 = 0.5*(E1+E3)
          else:
              print "Can't find satisfactory E in range"
              break
  
          i += 1
      x = inormalise_e(l,E2)
      if x[((l+1)**2)/0.005] &gt; (x[2*((l+1)**2)/0.005]) and iturns_e(l,E2+1e-20)==1:
          print 'Energy of state: n =',(l+1),'is: ',(E2*(10**9)),'eV'
          l += 1
      else:
          E1 = E2+10e-20
  </code></pre>