title: scipy.integrate.solve_ivp vectorized
body: <p>Trying to use the vectorized option for solve_ivp and strangely it throws an error that y0 must be 1 dimensional. 
  MWE :</p>
  
  <pre><code>from scipy.integrate import solve_ivp
  import numpy as np
  import math
  
  def f(t, y):
      theta = math.pi/4
      ham = np.array([[1,0],[1,np.exp(-1j*theta*t)]])
      return-1j * np.dot(ham,y)
  
  
  def main():
  
      y0 = np.eye(2,dtype= np.complex128)
      t0 = 0
      tmax = 10**(-6)
      sol=solve_ivp( lambda t,y :f(t,y),(t0,tmax),y0,method='RK45',vectorized=True)
      print(sol.y)
  
  if __name__ == '__main__':
      main()
  </code></pre>
  
  <blockquote>
    <p>The calling signature is fun(t, y). Here t is a scalar, and there are two options for the ndarray y: It can either have shape (n,); then fun must return array_like with shape (n,). Alternatively it can have shape (n, k); then fun must return an array_like with shape (n, k), i.e. each column corresponds to a single column in y. The choice between the two options is determined by vectorized argument (see below). The vectorized implementation allows a faster approximation of the Jacobian by finite differences (required for stiff solvers).</p>
  </blockquote>
  
  <p>Error :</p>
  
  <blockquote>
    <p>ValueError: <code>y0</code> must be 1-dimensional.</p>
  </blockquote>
  
  <p>Python 3.6.8</p>
  
  <p>scipy.<strong>version</strong> 
  '1.2.1'</p>