The `scipy.integrate.ode`

interface to integration routines provides a method for stopping the integration if a constraint is violated at any step, `set_solout`

. However, I cannot get this method to work, even in the simplest examples. Here’s one attempt:

```
import numpy as np
from scipy.integrate import ode
def f(t, y):
"""Exponential decay."""
return -y
def solout(t, y):
if y[0] < 0.5:
return -1
else:
return 0
y_initial = 1
t_initial = 0
r = ode(f).set_integrator('dopri5') # Integrator that supports solout
r.set_initial_value(y_initial, t_initial)
r.set_solout(solout)
# Integrate until t = 5, but stop when solout constraint violated
r.integrate(5)
# The time when solout should have terminated integration:
intersection_time = np.log(2)
```

The integration should have been stopped by solout when `t = log(2) = 0.693...`

, but instead happily continues until `t = 5`

, when `y = 0.007`

.

Is this a bug in `scipy`

, or am I not using `set_solout`

correctly?

Best answer

It turns out you need to call `set_solout`

*before* calling `set_initial_value`

. (I figured this out by studying the `set_solout`

tests in the `scipy`

test suite.) So, reversing the order of the two calls in my question code produces the correct result.

Even if this behavior is correct, it ought to be mentioned in the documentation for `set_solout`

. I’ve posted an issue with SciPy on GitHub.

**UPDATE:** This issue is fixed in SciPy 0.17.0; `set_solout`

will work even if called after `set_initial_value`

, and the question code will produce the correct result.