title: Restart python scipy.solve_ivp integrator with new parameters after terminal event occurs
body: <p>I'm using scipy.integrate.solve_ivp to solve a system of ODEs because it has the event functions. </p>
  
  <p>The reason why I need this function is that during the integration sometimes I get a singular matrix, and everytime that happens I need to finish the integration and restart it with new parameters. </p>
  
  <p>I would like to know if is possible to restart the scipy.integrate.solve_ivp with new parameters after a terminal event has occurred, and if so how could I do it. </p>
  
  <p>Any help would be very much appreciated.  </p>
  
  <p>This is my current script based on an example from 
  <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html" rel="nofollow noreferrer">https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html</a>:</p>
  
  <pre><code>#!/usr/bin/env python
  # -*- coding: utf-8 -*-
  
  import numpy as np
  from scipy.integrate import solve_ivp
  import matplotlib.pyplot as plt
  import matplotlib.animation as animation
  from animate_plot import animate
  
  def upward_cannon(t, y):
      return [y[1], -0.5]
  
  def hit_ground(t, y):
      return y[0]
  
  def apex(t,y):
      return y[1]
  
  hit_ground.terminal = True
  hit_ground.direction = -1
  
  t0 = 0
  
  tf = 10
  
  sol = solve_ivp(upward_cannon, t_span=[t0, tf], y0=[90, 10], t_eval=np.arange(t0,tf, 0.01), events=[hit_ground,apex],
                  dense_output=True)
  
  linesData = { 1: [[-0.0, 0.0],[0.0, 0.0]]}#,
              # 2: [[-0.5, 0],[0.5, 0.0]]}#, 3: [[-0.5, 0],[0.5, 0]]}
  
  pointsofInterest = {}#3: [[0.5, 0.0]]}#, 2: [[180.0, 10]]}
  
  model_markers = np.array([])
  
  plot_title = 'Upward Particle'
  plot_legend = ['Forward Dynamics']
  
  q_rep = sol.y.T[:,0]
  
  fig = plt.figure()
  
  ax = fig.add_subplot(111)
  
  xs = np.arange(t0, tf, 0.01)
  
  for idx in range(0,q_rep.shape[0]):   #looping statement;declare the total number of frames
      y=q_rep[idx]       # traveling Sine wave
      ax.cla()
      ax.scatter(xs[idx],y, s=50)
      plt.ylim([-10, 190])
      plt.xlim([-100, 100])
      plt.pause(0.001)
  plt.show()
  </code></pre>
  
  <p>Thank you in advance.</p>
  
  <p>Kind Regards</p>