must behave
    such that ``fun(t, y)[:, i] == fun(t, y[:, i])`` (i.e. each column of
    the returned array is the time derivative of the state corresponding
    with a column of ``y``).

    Setting ``vectorized=True`` allows for faster finite difference
    approximation of the Jacobian by methods 'Radau' and 'BDF', but
    will result in slower execution for this solver.

Attributes
----------
n : int
    Number of equations.
status : string
    Current status of the solver: 'running', 'finished' or 'failed'.
t_bound : float
    Boundary time.
direction : float
    Integration direction: +1 or -1.
t : float
    Current time.
y : ndarray
    Current state.
t_old : float
    Previous time. None if no steps were made yet.
nfev : int
    Number of evaluations of the right-hand side.
njev : int
    Number of evaluations of the Jacobian.

References
----------
.. [1] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
       Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
       pp. 55-64, 1983.
.. [2] L. Petzold, "Automatic selection of methods for solving stiff and
       nonstiff systems of ordinary differential equations", SIAM Journal
       on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
       1983.
Ng