title: Using solve_ivp instead of odeint to solve initial problem value
body: <p>Currently, I solve the following ODE system of equations using odeint</p>
  
  <p>dx/dt = (-x + u)/2.0</p>
  
  <p>dy/dt = (-y + x)/5.0</p>
  
  <p>initial conditions: x = 0, y = 0</p>
  
  <p>However, I would like to use solve_ivp which seems to be the recommended option for this type of problems, but honestly I don't know how to adapt the code...</p>
  
  <p>Here is the code I'm using with odeint:</p>
  
  <pre><code>import numpy as np
  from scipy.integrate import odeint, solve_ivp
  import matplotlib.pyplot as plt
  
  def model(z, t, u):
      x = z[0]
      y = z[1]
      dxdt = (-x + u)/2.0
      dydt = (-y + x)/5.0
      dzdt = [dxdt, dydt]
      return dzdt
  
  def main():
      # initial condition
      z0 = [0, 0]
  
      # number of time points
      n = 401
  
      # time points
      t = np.linspace(0, 40, n)
  
      # step input
      u = np.zeros(n)
      # change to 2.0 at time = 5.0
      u[51:] = 2.0
  
      # store solution
      x = np.empty_like(t)
      y = np.empty_like(t)
      # record initial conditions
      x[0] = z0[0]
      y[0] = z0[1]
  
      # solve ODE
      for i in range(1, n):
          # span for next time step
          tspan = [t[i-1], t[i]]
          # solve for next step
          z = odeint(model, z0, tspan, args=(u[i],))
          # store solution for plotting
          x[i] = z[1][0]
          y[i] = z[1][1]
          # next initial condition
          z0 = z[1]
  
      # plot results
      plt.plot(t,u,'g:',label='u(t)')
      plt.plot(t,x,'b-',label='x(t)')
      plt.plot(t,y,'r--',label='y(t)')
      plt.ylabel('values')
      plt.xlabel('time')
      plt.legend(loc='best')
      plt.show()
  
  main()
  </code></pre>