title: Solve_ivp integration gets stuck if initial condition triggers event with terminal = True
body: <p>I work in spacecraft trajectory design and I am implementing an algorithm for finding periodic orbits that needs to detect whenever a trajectory crosses a specific plane. I am using <code>scipy.integrate</code>'s <code>solve_ivp</code> and its <code>events</code> functionality. </p>
  
  <p>Ideally, I would like to count when events are triggered only stop the integration after a certain number of events are triggered, but since I couldn't find any similar option in <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html" rel="nofollow noreferrer">solve_ivp's documentation</a>, I am simply using a for cycle to chain these integration arcs and count the number of events. However, my problem is that <code>solve_ivp</code> gets stuck when the initial condition triggers the event and <code>event.terminal = True</code>. </p>
  
  <p>You can find an example code snipet below:</p>
  
  <pre><code>import numpy as np
  from scipy.integrate import solve_ivp
  
  def eq_motion(t,Y):
  
      Ydot = np.zeros((6,1))
      r = (Y[0,0]**2 + Y[1,0]**2 + Y[2,0]**2)**0.5
  
      Ydot[0] = Y[3,0]
      Ydot[1] = Y[4,0] 
      Ydot[2] = Y[5,0]
      Ydot[3] = 2*Y[4,0] + 3*Y[0,0] - Y[0,0]/r**3
      Ydot[4] = -2*Y[3,0] - Y[1,0]/r**3
      Ydot[5] = -Y[2,0] - Y[2,0]/r**3
  
      return Ydot
  
  def event_xz_plane(t,Y):
      return Y[1]
  
  event_xz_plane.terminal = True
  event_xz_plane.direction = 0
  
  # initial conditions
  Y0 = np.array([-0.006489114696252, 0, 0.107462057974196 + 0.0006, 0, 4.165257687202607,0])
  t_span = np.array([0, 1.892845635529040*2])
  
  n_plane_intersections = 4
  
  for i in np.arange(n_plane_intersections):
  
      integrate = solve_ivp(fun=eq_motion, t_span=t_span, y0=Y0, method = 'LSODA', 
                              vectorized = True, rtol = 1e-13, atol = 1e-14,
                              events = (event_xz_plane),dense_output=True)
  
      print('Number of integrations before event: {}'.format(integrate.y.shape[1]-1))
      print(integrate.y_events[0])
      print('\n')
  
      # attempting to chain the different integration arcs
      Y0 = integrate.y[:,-1]
  </code></pre>
  
  <p>The initial conditions chosen are close to those of a periodic orbit and produce the following <a href="https://i.stack.imgur.com/BDwii.png" rel="nofollow noreferrer">trajectory</a>, where you can see there are 4 intersections with the xz-plane, including the initial condition, which is on that plane.</p>
  
  <p>However, the above script prints out the following:</p>
  
  <pre><code>Number of integrations before event: 1
  [[-0.00648911  0.          0.10806206  0.          4.16525769  0.        ]]
  
  
  Number of integrations before event: 1
  [[-0.00648911  0.          0.10806206  0.          4.16525769  0.        ]]
  
  
  Number of integrations before event: 1
  [[-0.00648911  0.          0.10806206  0.          4.16525769  0.        ]]
  
  
  Number of integrations before event: 1
  [[-0.00648911  0.          0.10806206  0.          4.16525769  0.        ]]
  </code></pre>
  
  <p>where you can see that it basically gets stuck in the initial condition. This happens, of course, because the initial condition is situated at a terminal condition, however, there seems to be no option to "pass" the first iteration of the integration.  </p>
  
  <p>I tried writing the event function as </p>
  
  <pre><code>def event_xz_plane(t,Y):
      return Y[1] + 1e3*(t==0) 
  </code></pre>
  
  <p>so that it would ignore the initial state, but then it simply gets "stuck" on the next plane intersection, and either way it doesn't seem like a good solution.</p>
  
  <p>The only other solution that I thought of was to set <code>terminate=False</code> on the event function and check the number of events triggered and the state vector at these instances. However, this is not ideal for my implementation since it depends on the <code>t_span</code> time set, which is an unknown in the algorithm. That is, if I need to get and <code>n</code> number of plane intersections, I would have to iteratively adjust t_span until I reach that number, instead of just allowing the integration to run until <code>n</code> plane intersections are triggered.</p>
  
  <p>Is there anyway to overcome this? </p>