title: scipy solve_ivp with negative first_step
body: <p>I want to solve a system of ODEs with scipy's solve_ivp and I need a negative first_step in order to evolve the solution, but </p>
  
  <pre><code>Nsol = solve_ivp(derivs, (N , Nend), ydoub, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, first_step=-1e-8)
  </code></pre>
  
  <p>returns</p>
  
  <pre><code>ValueError: `first_step` must be positive.
  </code></pre>
  
  <p>Any ideas how to resolve this or find a way around it?</p>
  
  <p>EDIT: here is the code which yields this:</p>
  
  <pre><code>import numpy
  from scipy.integrate import solve_ivp
  
  w = numpy.array ( [  0.00000000e+00,   1.00000000e+00,   3.17214587e-01,
          -3.41988549e-01,  -1.50137165e-05,  -2.48117074e-02,
           1.17624224e-03,  -1.27149037e-04] )
  
  def derivs2 (t, w):
      dydN = numpy.zeros(2 , dtype=float , order='C')
  
      dydN[0] = 0.0
  
      dydN[1] = y[1] * y[2]
      dydN[2] = y[2] * ( y[3] + 2.0 * y[2] )
      dydN[3] = 2.0 * y[4] - 5.0 * y[2] * y[3] - 12.0 * y[2] * y[2]
  
      for i in range (4 , NEQS-1):
          dydN[i] = ( 0.5 * (i-3) * y[3] + (i-4) * y[2] ) * y[i] + y[i+1]
  
      dydN[NEQS-1] = ( 0.5 * (NEQS-4) * y[3] + (NEQS-5) * y[2] ) * y[NEQS-1]
  
      return dydN
  
  Nsol = solve_ivp(derivs, (1000.0 , 0.0), w, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False , first_step=-1e-6)
  </code></pre>