title: Debugging, solve_ivp when used with toms748
body: <p>I am facing an error when implementing the following code</p>
  <pre><code>from scipy.optimize import toms748
  from scipy.integrate import solve_ivp
  
  def f(r):
      return lambda x: x-r
  
  def E(t,r):
      return -toms748(f(r),r-1,r+1)
  
  sol=solve_ivp(E,(0,10),[1])
  </code></pre>
  <p>The error obtained is as follows</p>
  <pre><code>Traceback (most recent call last):
  File &quot;C:\Users\User\Documents\Project codes\Density.py&quot;, line 10, in &lt;module&gt;
  sol=solve_ivp(E,(0,10),[1])
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\integrate\_ivp\ivp.py&quot;, line 576, in solve_ivp
  message = solver.step()
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\integrate\_ivp\base.py&quot;, line 181, in step
  success, message = self._step_impl()
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\integrate\_ivp\rk.py&quot;, line 144, in _step_impl
  y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A,
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\integrate\_ivp\rk.py&quot;, line 64, in rk_step
  K[s] = fun(t + c * h, y + dy)
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\integrate\_ivp\base.py&quot;, line 138, in fun
  return self.fun_single(t, y)
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\integrate\_ivp\base.py&quot;, line 20, in fun_wrapped
  return np.asarray(fun(t, y), dtype=dtype)
  File &quot;C:\Users\User\Documents\Project codes\Density.py&quot;, line 8, in E
  return -toms748(f(r),r-1,r+1)
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\optimize\zeros.py&quot;, line 1361, in toms748
  result = solver.solve(f, a, b, args=args, k=k, xtol=xtol, rtol=rtol,
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\optimize\zeros.py&quot;, line 1225, in solve
  status, xn = self.iterate()
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\optimize\zeros.py&quot;, line 1144, in iterate
  c = _newton_quadratic(self.ab, self.fab, d, fd, nsteps)
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\optimize\zeros.py&quot;, line 1004, in _newton_quadratic
  _, B, A = _compute_divided_differences([a, b, d], [fa, fb, fd],
  File &quot;C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
  packages\scipy\optimize\zeros.py&quot;, line 959, in _compute_divided_differences
  row = np.diff(row)[:] / denom
  ValueError: operands could not be broadcast together with shapes (3,0) (2,1)
  </code></pre>
  <p>toms748 is a root finding alogrithm that takes in a callable function and two scalars that bound the values between which the root is searched for.
  Thus E(t,r) is just E(t,r)=-r and the differential equation implemented above is dr/dt=-r with initial condition r(0)=1. The solution is just r(t)=exp(-t).</p>
  <p>Now the thing that perplexes me even more is, when I remove the minus sign from E(t,r) i.e let</p>
  <pre><code>def E(t,r):
      return toms748(f(r),r-1,r+1)
  </code></pre>
  <p>then the code runs just fine and no errors are thrown up.</p>
  <p>The above are all toy functions. The actual implementation is much more complicated. The above is the simplest code I could obtain that gives the same error.</p>