title: How can I control odeint to stop integration when the result reach a threshold?
body: <p>Here is my code.</p>
  
  <pre><code>import numpy as np
  from scipy.integrate import odeint
  
  #Constant
  R0=1.475
  gamma=2.
  ScaleMeVfm3toEskm3 = 8.92*np.power(10.,-7.)
  
  def EOSe(p):
      return np.power((p/450.785),(1./gamma))
  
  def M(m,r):
      return (4./3.)*np.pi*np.power(r,3.)*p
  
  # function that returns dz/dt
  def model(z,r):
      p, m = z
      dpdr = -((R0*EOSe(p)*m)/(np.power(r,2.)))*(1+(p/EOSe(p)))*(1+((4*math.pi*(np.power(r,3))*p)/(m)))*((1-((2*R0)*m)/(r))**(-1.))
      dmdr = 4.*math.pi*(r**2.)*EOSe(p)
      dzdr = [dpdr,dmdr]
      return dzdr
  
  # initial condition
  r0=10.**-12.
  p0=10**-6.
  z0 = [p0, M(r0, p0)]
  
  # radius
  r = np.linspace(r0, 15, 100000)
  
  # solve ODE
  z = odeint(model,z0,r)
  </code></pre>
  
  <p>The result of <code>z[:,0]</code> keeps decreasing as I expected. But what I want is only positive values. One may run the code and try <code>print(z[69306])</code> and it will show <code>[2.89636405e-11 5.46983202e-01]</code>. That is the last point I want the <code>odeint</code> to stop integration.</p>
  
  <p>Of course, the provided code shows</p>
  
  <pre><code>RuntimeWarning: invalid value encountered in power
    return np.power((p/450.785),(1./gamma))
  </code></pre>
  
  <p>because the result of p starts being negative. For any further points, the <code>odeint</code> yields the result <code>[nan nan]</code>. </p>
  
  <p>However, I could use np.nanmin() to find the minimum of z[:,0] that is not <code>nan</code>. But I have a set of <code>p0</code> values for my work. I will need to call <code>odeint</code> in a loop like</p>
  
  <pre><code>P=np.linspace(10**-8.,10**-2.,10000)
  for p0 in P:
  #the code for solving ode provided above.
  </code></pre>
  
  <p>which takes more time.</p>
  
  <p>I think it would reduce a time for execution if I can just stop at before z[:,0] going to be negative a value?</p>