title: Solving ODE with SciPy gives invalid value encountered in double_scalars error
body: <p>I'm trying to solve the Lane-Emden equation for arbitrary values of n (polytropic index). In order to use SciPy, I have represented the second-order ODE as a set of two coupled first-order ODEs. I have the following code:</p>
  
  <pre><code>import matplotlib.pyplot as plt
  import numpy as np
  from scipy.integrate import odeint
  
  def poly(y,xi,n):  
      theta, phi = y  
      dydt = [phi/(xi**2), -(xi**2)*theta**n]  
      return dydt
  </code></pre>
  
  <p>Here, I defined phi = theta' </p>
  
  <pre><code>y0 = [1.,0.]  
  xi = np.linspace(0., 16., 201)    
  for n in range(0,11):  
      sol = odeint(poly, y0, xi, args=(n/2.,))  
      plt.plot(xi, sol[:, 0], label=str(n/2.))  
  plt.legend(loc='best')  
  plt.xlabel('xi')  
  plt.grid()  
  plt.show()
  </code></pre>
  
  <p>However, running this code results in the following error:</p>
  
  <blockquote>
    <p>RuntimeWarning: invalid value encountered in double_scalars<br>
    app.launch_new_instance()</p>
  </blockquote>
  
  <p>At first, I thought this was a result of the singularity of the equation at xi = 0, so I changed my integration interval:</p>
  
  <pre><code>xi = np.linspace(1e-10, 16., 201)
  </code></pre>
  
  <p>This fixes the problem at n = 0., but not at other values of n. The plots I get don't make any sense and are just plain wrong.</p>
  
  <p>Why do I get this error message, and how can I fix my code?</p>