We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
lsoda fails to detect stiff problem when called from solve_ivp
<!-- 
Thank you for taking the time to report a SciPy issue.
Please describe the issue in detail, and for bug reports
fill in the fields below. You can delete the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
-->

Both the `ode` and the `solve_ivp` solver interfaces wrap the venerable lsoda solver.  As far as I can see, `solve_ivp` even internally calls the `ode` wrapper.  Strangely enough, when tested on a stiff problem, lsoda when called through `ode` handles stiffness just fine.  However, when called through `solve_ivp`, it behaves as if using the non-stiff solver.  The only parameter difference I can spot is the reduced default `atol` and `rtol` parameter, however, even setting them to the `ode` defaults will not change the behavior of `solve_ivp`.  (Off-topic comment: I find the default tolerances for `solve_ivp` unnecessarily poor, especially given that any other blackbox solver I know, including the two older Scipy APIs, have stricter defaults.  I would strongly recommend using common defaults across Scipy, otherwise an already confusing situation is made more confusing and less robust.)

### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be marked as a code block automatically 
-->


```python
#!/usr/bin/python3

from pylab import *
from scipy.integrate import solve_ivp
from scipy.integrate import ode
import time
    
mu = 1000
y0 = r_[0.5, 0.5]
T = 500
tt = linspace(T/200,T,200)

def van_der_pol(t, y):
    return r_[y[1], mu*(1.0-y[0]**2)*y[1]-y[0]]

c1 = time.process_time()
sol1 = solve_ivp(van_der_pol, [0,T], y0, t_eval=tt,
                 method = 'LSODA', rtol=1e-6, atol=1e-12)

c2 = time.process_time()
r2 = ode(van_der_pol).set_integrator('lsoda')
r2.set_initial_value(y0)
sol2 = array([r2.integrate(T) for T in tt])

c3 = time.process_time() 
print ("Time for LSODA (solve_ivp): " + str(c2-c1))
print ("Time for lsoda (ode):       " + str(c3-c2))

figure()
plot(sol1.t, sol1.y[0,:], label='LSODA (solve_ivp)')
plot(tt, sol2[:,0], label='lsoda (ode)')
legend()
show()
```
This yields the following output (graphical output is qualitatively OK with both cases being very close):
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
Time for LSODA (solve_ivp): 35.811249636000014
Time for lsoda (ode):       0.024488598000004913
```
So there is more than a factor 1000 performance difference, these numbers indicates that in the ill-performing case, `lsoda` has not switched to using a BDF method.

I tried to find any substantial differences in the code paths of the two wrappers, without success.  

### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
-->
1.1.0 1.15.1 sys.version_info(major=3, minor=7, micro=2, releaselevel='final', serial=0)


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `scipy/integrate/_ivp/lsoda.py`
Content:
```
1 import numpy as np
2 from scipy.integrate import ode
3 from .common import validate_tol, validate_first_step, warn_extraneous
4 from .base import OdeSolver, DenseOutput
5 
6 
7 class LSODA(OdeSolver):
8     """Adams/BDF method with automatic stiffness detection and switching.
9 
10     This is a wrapper to the Fortran solver from ODEPACK [1]_. It switches
11     automatically between the nonstiff Adams method and the stiff BDF method.
12     The method was originally detailed in [2]_.
13 
14     Parameters
15     ----------
16     fun : callable
17         Right-hand side of the system. The calling signature is ``fun(t, y)``.
18         Here ``t`` is a scalar, and there are two options for the ndarray ``y``:
19         It can either have shape (n,); then ``fun`` must return array_like with
20         shape (n,). Alternatively it can have shape (n, k); then ``fun``
21         must return an array_like with shape (n, k), i.e. each column
22         corresponds to a single column in ``y``. The choice between the two
23         options is determined by `vectorized` argument (see below). The
24         vectorized implementation allows a faster approximation of the Jacobian
25         by finite differences (required for this solver).
26     t0 : float
27         Initial time.
28     y0 : array_like, shape (n,)
29         Initial state.
30     t_bound : float
31         Boundary time - the integration won't continue beyond it. It also
32         determines the direction of the integration.
33     first_step : float or None, optional
34         Initial step size. Default is ``None`` which means that the algorithm
35         should choose.
36     min_step : float, optional
37         Minimum allowed step size. Default is 0.0, i.e. the step size is not
38         bounded and determined solely by the solver.
39     max_step : float, optional
40         Maximum allowed step size. Default is np.inf, i.e. the step size is not
41         bounded and determined solely by the solver.
42     rtol, atol : float and array_like, optional
43         Relative and absolute tolerances. The solver keeps the local error
44         estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
45         relative accuracy (number of correct digits). But if a component of `y`
46         is approximately below `atol`, the error only needs to fall within
47         the same `atol` threshold, and the number of correct digits is not
48         guaranteed. If components of y have different scales, it might be
49         beneficial to set different `atol` values for different components by
50         passing array_like with shape (n,) for `atol`. Default values are
51         1e-3 for `rtol` and 1e-6 for `atol`.
52     jac : None or callable, optional
53         Jacobian matrix of the right-hand side of the system with respect to
54         ``y``. The Jacobian matrix has shape (n, n) and its element (i, j) is
55         equal to ``d f_i / d y_j``. The function will be called as
56         ``jac(t, y)``. If None (default), the Jacobian will be
57         approximated by finite differences. It is generally recommended to
58         provide the Jacobian rather than relying on a finite-difference
59         approximation.
60     lband, uband : int or None
61         Parameters defining the bandwidth of the Jacobian,
62         i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``. Setting
63         these requires your jac routine to return the Jacobian in the packed format:
64         the returned array must have ``n`` columns and ``uband + lband + 1``
65         rows in which Jacobian diagonals are written. Specifically
66         ``jac_packed[uband + i - j , j] = jac[i, j]``. The same format is used
67         in `scipy.linalg.solve_banded` (check for an illustration).
68         These parameters can be also used with ``jac=None`` to reduce the
69         number of Jacobian elements estimated by finite differences.
70     vectorized : bool, optional
71         Whether `fun` is implemented in a vectorized fashion. A vectorized
72         implementation offers no advantages for this solver. Default is False.
73 
74     Attributes
75     ----------
76     n : int
77         Number of equations.
78     status : string
79         Current status of the solver: 'running', 'finished' or 'failed'.
80     t_bound : float
81         Boundary time.
82     direction : float
83         Integration direction: +1 or -1.
84     t : float
85         Current time.
86     y : ndarray
87         Current state.
88     t_old : float
89         Previous time. None if no steps were made yet.
90     nfev : int
91         Number of evaluations of the right-hand side.
92     njev : int
93         Number of evaluations of the Jacobian.
94 
95     References
96     ----------
97     .. [1] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
98            Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
99            pp. 55-64, 1983.
100     .. [2] L. Petzold, "Automatic selection of methods for solving stiff and
101            nonstiff systems of ordinary differential equations", SIAM Journal
102            on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
103            1983.
104     """
105     def __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0,
106                  max_step=np.inf, rtol=1e-3, atol=1e-6, jac=None, lband=None,
107                  uband=None, vectorized=False, **extraneous):
108         warn_extraneous(extraneous)
109         super(LSODA, self).__init__(fun, t0, y0, t_bound, vectorized)
110 
111         if first_step is None:
112             first_step = 0  # LSODA value for automatic selection.
113         else:
114             first_step = validate_first_step(first_step, t0, t_bound)
115 
116         first_step *= self.direction
117 
118         if max_step == np.inf:
119             max_step = 0  # LSODA value for infinity.
120         elif max_step <= 0:
121             raise ValueError("`max_step` must be positive.")
122 
123         if min_step < 0:
124             raise ValueError("`min_step` must be nonnegative.")
125 
126         rtol, atol = validate_tol(rtol, atol, self.n)
127 
128         if jac is None:  # No lambda as PEP8 insists.
129             def jac():
130                 return None
131 
132         solver = ode(self.fun, jac)
133         solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,
134                               min_step=min_step, first_step=first_step,
135                               lband=lband, uband=uband)
136         solver.set_initial_value(y0, t0)
137 
138         # Inject t_bound into rwork array as needed for itask=5.
139         solver._integrator.rwork[0] = self.t_bound
140         solver._integrator.call_args[4] = solver._integrator.rwork
141 
142         self._lsoda_solver = solver
143 
144     def _step_impl(self):
145         solver = self._lsoda_solver
146         integrator = solver._integrator
147 
148         # From lsoda.step and lsoda.integrate itask=5 means take a single
149         # step and do not go past t_bound.
150         itask = integrator.call_args[2]
151         integrator.call_args[2] = 5
152         solver._y, solver.t = integrator.run(
153             solver.f, solver.jac, solver._y, solver.t,
154             self.t_bound, solver.f_params, solver.jac_params)
155         integrator.call_args[2] = itask
156 
157         if solver.successful():
158             self.t = solver.t
159             self.y = solver._y
160             # From LSODA Fortran source njev is equal to nlu.
161             self.njev = integrator.iwork[12]
162             self.nlu = integrator.iwork[12]
163             return True, None
164         else:
165             return False, 'Unexpected istate in LSODA.'
166 
167     def _dense_output_impl(self):
168         iwork = self._lsoda_solver._integrator.iwork
169         rwork = self._lsoda_solver._integrator.rwork
170 
171         order = iwork[14]
172         h = rwork[11]
173         yh = np.reshape(rwork[20:20 + (order + 1) * self.n],
174                         (self.n, order + 1), order='F').copy()
175 
176         return LsodaDenseOutput(self.t_old, self.t, h, order, yh)
177 
178 
179 class LsodaDenseOutput(DenseOutput):
180     def __init__(self, t_old, t, h, order, yh):
181         super(LsodaDenseOutput, self).__init__(t_old, t)
182         self.h = h
183         self.yh = yh
184         self.p = np.arange(order + 1)
185 
186     def _call_impl(self, t):
187         if t.ndim == 0:
188             x = ((t - self.t) / self.h) ** self.p
189         else:
190             x = ((t - self.t) / self.h) ** self.p[:, None]
191 
192         return np.dot(self.yh, x)
193
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/scipy/integrate/_ivp/lsoda.py b/scipy/integrate/_ivp/lsoda.py
--- a/scipy/integrate/_ivp/lsoda.py
+++ b/scipy/integrate/_ivp/lsoda.py
@@ -125,10 +125,6 @@
 
         rtol, atol = validate_tol(rtol, atol, self.n)
 
-        if jac is None:  # No lambda as PEP8 insists.
-            def jac():
-                return None
-
         solver = ode(self.fun, jac)
         solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,
                               min_step=min_step, first_step=first_step,
@@ -150,7 +146,7 @@
         itask = integrator.call_args[2]
         integrator.call_args[2] = 5
         solver._y, solver.t = integrator.run(
-            solver.f, solver.jac, solver._y, solver.t,
+            solver.f, solver.jac or (lambda: None), solver._y, solver.t,
             self.t_bound, solver.f_params, solver.jac_params)
         integrator.call_args[2] = itask

{"golden_diff": "diff --git a/scipy/integrate/_ivp/lsoda.py b/scipy/integrate/_ivp/lsoda.py\n--- a/scipy/integrate/_ivp/lsoda.py\n+++ b/scipy/integrate/_ivp/lsoda.py\n@@ -125,10 +125,6 @@\n \n         rtol, atol = validate_tol(rtol, atol, self.n)\n \n-        if jac is None:  # No lambda as PEP8 insists.\n-            def jac():\n-                return None\n-\n         solver = ode(self.fun, jac)\n         solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,\n                               min_step=min_step, first_step=first_step,\n@@ -150,7 +146,7 @@\n         itask = integrator.call_args[2]\n         integrator.call_args[2] = 5\n         solver._y, solver.t = integrator.run(\n-            solver.f, solver.jac, solver._y, solver.t,\n+            solver.f, solver.jac or (lambda: None), solver._y, solver.t,\n             self.t_bound, solver.f_params, solver.jac_params)\n         integrator.call_args[2] = itask\n", "issue": "lsoda fails to detect stiff problem when called from solve_ivp\n<!-- \r\nThank you for taking the time to report a SciPy issue.\r\nPlease describe the issue in detail, and for bug reports\r\nfill in the fields below. You can delete the sections that \r\ndon't apply to your issue. You can view the final output\r\nby clicking the preview button above.\r\n-->\r\n\r\nBoth the `ode` and the `solve_ivp` solver interfaces wrap the venerable lsoda solver.  As far as I can see, `solve_ivp` even internally calls the `ode` wrapper.  Strangely enough, when tested on a stiff problem, lsoda when called through `ode` handles stiffness just fine.  However, when called through `solve_ivp`, it behaves as if using the non-stiff solver.  The only parameter difference I can spot is the reduced default `atol` and `rtol` parameter, however, even setting them to the `ode` defaults will not change the behavior of `solve_ivp`.  (Off-topic comment: I find the default tolerances for `solve_ivp` unnecessarily poor, especially given that any other blackbox solver I know, including the two older Scipy APIs, have stricter defaults.  I would strongly recommend using common defaults across Scipy, otherwise an already confusing situation is made more confusing and less robust.)\r\n\r\n### Reproducing code example:\r\n<!-- \r\nIf you place your code between the triple backticks below, \r\nit will be marked as a code block automatically \r\n-->\r\n\r\n\r\n```python\r\n#!/usr/bin/python3\r\n\r\nfrom pylab import *\r\nfrom scipy.integrate import solve_ivp\r\nfrom scipy.integrate import ode\r\nimport time\r\n    \r\nmu = 1000\r\ny0 = r_[0.5, 0.5]\r\nT = 500\r\ntt = linspace(T/200,T,200)\r\n\r\ndef van_der_pol(t, y):\r\n    return r_[y[1], mu*(1.0-y[0]**2)*y[1]-y[0]]\r\n\r\nc1 = time.process_time()\r\nsol1 = solve_ivp(van_der_pol, [0,T], y0, t_eval=tt,\r\n                 method = 'LSODA', rtol=1e-6, atol=1e-12)\r\n\r\nc2 = time.process_time()\r\nr2 = ode(van_der_pol).set_integrator('lsoda')\r\nr2.set_initial_value(y0)\r\nsol2 = array([r2.integrate(T) for T in tt])\r\n\r\nc3 = time.process_time() \r\nprint (\"Time for LSODA (solve_ivp): \" + str(c2-c1))\r\nprint (\"Time for lsoda (ode):       \" + str(c3-c2))\r\n\r\nfigure()\r\nplot(sol1.t, sol1.y[0,:], label='LSODA (solve_ivp)')\r\nplot(tt, sol2[:,0], label='lsoda (ode)')\r\nlegend()\r\nshow()\r\n```\r\nThis yields the following output (graphical output is qualitatively OK with both cases being very close):\r\n<!-- If any, paste the *full* error message inside a code block\r\nas above (starting from line Traceback)\r\n-->\r\n\r\n```\r\nTime for LSODA (solve_ivp): 35.811249636000014\r\nTime for lsoda (ode):       0.024488598000004913\r\n```\r\nSo there is more than a factor 1000 performance difference, these numbers indicates that in the ill-performing case, `lsoda` has not switched to using a BDF method.\r\n\r\nI tried to find any substantial differences in the code paths of the two wrappers, without success.  \r\n\r\n### Scipy/Numpy/Python version information:\r\n<!-- You can simply run the following and paste the result in a code block\r\n```\r\nimport sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)\r\n```\r\n-->\r\n1.1.0 1.15.1 sys.version_info(major=3, minor=7, micro=2, releaselevel='final', serial=0)\r\n\n", "before_files": [{"content": "import numpy as np\nfrom scipy.integrate import ode\nfrom .common import validate_tol, validate_first_step, warn_extraneous\nfrom .base import OdeSolver, DenseOutput\n\n\nclass LSODA(OdeSolver):\n    \"\"\"Adams/BDF method with automatic stiffness detection and switching.\n\n    This is a wrapper to the Fortran solver from ODEPACK [1]_. It switches\n    automatically between the nonstiff Adams method and the stiff BDF method.\n    The method was originally detailed in [2]_.\n\n    Parameters\n    ----------\n    fun : callable\n        Right-hand side of the system. The calling signature is ``fun(t, y)``.\n        Here ``t`` is a scalar, and there are two options for the ndarray ``y``:\n        It can either have shape (n,); then ``fun`` must return array_like with\n        shape (n,). Alternatively it can have shape (n, k); then ``fun``\n        must return an array_like with shape (n, k), i.e. each column\n        corresponds to a single column in ``y``. The choice between the two\n        options is determined by `vectorized` argument (see below). The\n        vectorized implementation allows a faster approximation of the Jacobian\n        by finite differences (required for this solver).\n    t0 : float\n        Initial time.\n    y0 : array_like, shape (n,)\n        Initial state.\n    t_bound : float\n        Boundary time - the integration won't continue beyond it. It also\n        determines the direction of the integration.\n    first_step : float or None, optional\n        Initial step size. Default is ``None`` which means that the algorithm\n        should choose.\n    min_step : float, optional\n        Minimum allowed step size. Default is 0.0, i.e. the step size is not\n        bounded and determined solely by the solver.\n    max_step : float, optional\n        Maximum allowed step size. Default is np.inf, i.e. the step size is not\n        bounded and determined solely by the solver.\n    rtol, atol : float and array_like, optional\n        Relative and absolute tolerances. The solver keeps the local error\n        estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a\n        relative accuracy (number of correct digits). But if a component of `y`\n        is approximately below `atol`, the error only needs to fall within\n        the same `atol` threshold, and the number of correct digits is not\n        guaranteed. If components of y have different scales, it might be\n        beneficial to set different `atol` values for different components by\n        passing array_like with shape (n,) for `atol`. Default values are\n        1e-3 for `rtol` and 1e-6 for `atol`.\n    jac : None or callable, optional\n        Jacobian matrix of the right-hand side of the system with respect to\n        ``y``. The Jacobian matrix has shape (n, n) and its element (i, j) is\n        equal to ``d f_i / d y_j``. The function will be called as\n        ``jac(t, y)``. If None (default), the Jacobian will be\n        approximated by finite differences. It is generally recommended to\n        provide the Jacobian rather than relying on a finite-difference\n        approximation.\n    lband, uband : int or None\n        Parameters defining the bandwidth of the Jacobian,\n        i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``. Setting\n        these requires your jac routine to return the Jacobian in the packed format:\n        the returned array must have ``n`` columns and ``uband + lband + 1``\n        rows in which Jacobian diagonals are written. Specifically\n        ``jac_packed[uband + i - j , j] = jac[i, j]``. The same format is used\n        in `scipy.linalg.solve_banded` (check for an illustration).\n        These parameters can be also used with ``jac=None`` to reduce the\n        number of Jacobian elements estimated by finite differences.\n    vectorized : bool, optional\n        Whether `fun` is implemented in a vectorized fashion. A vectorized\n        implementation offers no advantages for this solver. Default is False.\n\n    Attributes\n    ----------\n    n : int\n        Number of equations.\n    status : string\n        Current status of the solver: 'running', 'finished' or 'failed'.\n    t_bound : float\n        Boundary time.\n    direction : float\n        Integration direction: +1 or -1.\n    t : float\n        Current time.\n    y : ndarray\n        Current state.\n    t_old : float\n        Previous time. None if no steps were made yet.\n    nfev : int\n        Number of evaluations of the right-hand side.\n    njev : int\n        Number of evaluations of the Jacobian.\n\n    References\n    ----------\n    .. [1] A. C. Hindmarsh, \"ODEPACK, A Systematized Collection of ODE\n           Solvers,\" IMACS Transactions on Scientific Computation, Vol 1.,\n           pp. 55-64, 1983.\n    .. [2] L. Petzold, \"Automatic selection of methods for solving stiff and\n           nonstiff systems of ordinary differential equations\", SIAM Journal\n           on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,\n           1983.\n    \"\"\"\n    def __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0,\n                 max_step=np.inf, rtol=1e-3, atol=1e-6, jac=None, lband=None,\n                 uband=None, vectorized=False, **extraneous):\n        warn_extraneous(extraneous)\n        super(LSODA, self).__init__(fun, t0, y0, t_bound, vectorized)\n\n        if first_step is None:\n            first_step = 0  # LSODA value for automatic selection.\n        else:\n            first_step = validate_first_step(first_step, t0, t_bound)\n\n        first_step *= self.direction\n\n        if max_step == np.inf:\n            max_step = 0  # LSODA value for infinity.\n        elif max_step <= 0:\n            raise ValueError(\"`max_step` must be positive.\")\n\n        if min_step < 0:\n            raise ValueError(\"`min_step` must be nonnegative.\")\n\n        rtol, atol = validate_tol(rtol, atol, self.n)\n\n        if jac is None:  # No lambda as PEP8 insists.\n            def jac():\n                return None\n\n        solver = ode(self.fun, jac)\n        solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,\n                              min_step=min_step, first_step=first_step,\n                              lband=lband, uband=uband)\n        solver.set_initial_value(y0, t0)\n\n        # Inject t_bound into rwork array as needed for itask=5.\n        solver._integrator.rwork[0] = self.t_bound\n        solver._integrator.call_args[4] = solver._integrator.rwork\n\n        self._lsoda_solver = solver\n\n    def _step_impl(self):\n        solver = self._lsoda_solver\n        integrator = solver._integrator\n\n        # From lsoda.step and lsoda.integrate itask=5 means take a single\n        # step and do not go past t_bound.\n        itask = integrator.call_args[2]\n        integrator.call_args[2] = 5\n        solver._y, solver.t = integrator.run(\n            solver.f, solver.jac, solver._y, solver.t,\n            self.t_bound, solver.f_params, solver.jac_params)\n        integrator.call_args[2] = itask\n\n        if solver.successful():\n            self.t = solver.t\n            self.y = solver._y\n            # From LSODA Fortran source njev is equal to nlu.\n            self.njev = integrator.iwork[12]\n            self.nlu = integrator.iwork[12]\n            return True, None\n        else:\n            return False, 'Unexpected istate in LSODA.'\n\n    def _dense_output_impl(self):\n        iwork = self._lsoda_solver._integrator.iwork\n        rwork = self._lsoda_solver._integrator.rwork\n\n        order = iwork[14]\n        h = rwork[11]\n        yh = np.reshape(rwork[20:20 + (order + 1) * self.n],\n                        (self.n, order + 1), order='F').copy()\n\n        return LsodaDenseOutput(self.t_old, self.t, h, order, yh)\n\n\nclass LsodaDenseOutput(DenseOutput):\n    def __init__(self, t_old, t, h, order, yh):\n        super(LsodaDenseOutput, self).__init__(t_old, t)\n        self.h = h\n        self.yh = yh\n        self.p = np.arange(order + 1)\n\n    def _call_impl(self, t):\n        if t.ndim == 0:\n            x = ((t - self.t) / self.h) ** self.p\n        else:\n            x = ((t - self.t) / self.h) ** self.p[:, None]\n\n        return np.dot(self.yh, x)\n", "path": "scipy/integrate/_ivp/lsoda.py"}], "after_files": [{"content": "import numpy as np\nfrom scipy.integrate import ode\nfrom .common import validate_tol, validate_first_step, warn_extraneous\nfrom .base import OdeSolver, DenseOutput\n\n\nclass LSODA(OdeSolver):\n    \"\"\"Adams/BDF method with automatic stiffness detection and switching.\n\n    This is a wrapper to the Fortran solver from ODEPACK [1]_. It switches\n    automatically between the nonstiff Adams method and the stiff BDF method.\n    The method was originally detailed in [2]_.\n\n    Parameters\n    ----------\n    fun : callable\n        Right-hand side of the system. The calling signature is ``fun(t, y)``.\n        Here ``t`` is a scalar, and there are two options for the ndarray ``y``:\n        It can either have shape (n,); then ``fun`` must return array_like with\n        shape (n,). Alternatively it can have shape (n, k); then ``fun``\n        must return an array_like with shape (n, k), i.e. each column\n        corresponds to a single column in ``y``. The choice between the two\n        options is determined by `vectorized` argument (see below). The\n        vectorized implementation allows a faster approximation of the Jacobian\n        by finite differences (required for this solver).\n    t0 : float\n        Initial time.\n    y0 : array_like, shape (n,)\n        Initial state.\n    t_bound : float\n        Boundary time - the integration won't continue beyond it. It also\n        determines the direction of the integration.\n    first_step : float or None, optional\n        Initial step size. Default is ``None`` which means that the algorithm\n        should choose.\n    min_step : float, optional\n        Minimum allowed step size. Default is 0.0, i.e. the step size is not\n        bounded and determined solely by the solver.\n    max_step : float, optional\n        Maximum allowed step size. Default is np.inf, i.e. the step size is not\n        bounded and determined solely by the solver.\n    rtol, atol : float and array_like, optional\n        Relative and absolute tolerances. The solver keeps the local error\n        estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a\n        relative accuracy (number of correct digits). But if a component of `y`\n        is approximately below `atol`, the error only needs to fall within\n        the same `atol` threshold, and the number of correct digits is not\n        guaranteed. If components of y have different scales, it might be\n        beneficial to set different `atol` values for different components by\n        passing array_like with shape (n,) for `atol`. Default values are\n        1e-3 for `rtol` and 1e-6 for `atol`.\n    jac : None or callable, optional\n        Jacobian matrix of the right-hand side of the system with respect to\n        ``y``. The Jacobian matrix has shape (n, n) and its element (i, j) is\n        equal to ``d f_i / d y_j``. The function will be called as\n        ``jac(t, y)``. If None (default), the Jacobian will be\n        approximated by finite differences. It is generally recommended to\n        provide the Jacobian rather than relying on a finite-difference\n        approximation.\n    lband, uband : int or None\n        Parameters defining the bandwidth of the Jacobian,\n        i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``. Setting\n        these requires your jac routine to return the Jacobian in the packed format:\n        the returned array must have ``n`` columns and ``uband + lband + 1``\n        rows in which Jacobian diagonals are written. Specifically\n        ``jac_packed[uband + i - j , j] = jac[i, j]``. The same format is used\n        in `scipy.linalg.solve_banded` (check for an illustration).\n        These parameters can be also used with ``jac=None`` to reduce the\n        number of Jacobian elements estimated by finite differences.\n    vectorized : bool, optional\n        Whether `fun` is implemented in a vectorized fashion. A vectorized\n        implementation offers no advantages for this solver. Default is False.\n\n    Attributes\n    ----------\n    n : int\n        Number of equations.\n    status : string\n        Current status of the solver: 'running', 'finished' or 'failed'.\n    t_bound : float\n        Boundary time.\n    direction : float\n        Integration direction: +1 or -1.\n    t : float\n        Current time.\n    y : ndarray\n        Current state.\n    t_old : float\n        Previous time. None if no steps were made yet.\n    nfev : int\n        Number of evaluations of the right-hand side.\n    njev : int\n        Number of evaluations of the Jacobian.\n\n    References\n    ----------\n    .. [1] A. C. Hindmarsh, \"ODEPACK, A Systematized Collection of ODE\n           Solvers,\" IMACS Transactions on Scientific Computation, Vol 1.,\n           pp. 55-64, 1983.\n    .. [2] L. Petzold, \"Automatic selection of methods for solving stiff and\n           nonstiff systems of ordinary differential equations\", SIAM Journal\n           on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,\n           1983.\n    \"\"\"\n    def __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0,\n                 max_step=np.inf, rtol=1e-3, atol=1e-6, jac=None, lband=None,\n                 uband=None, vectorized=False, **extraneous):\n        warn_extraneous(extraneous)\n        super(LSODA, self).__init__(fun, t0, y0, t_bound, vectorized)\n\n        if first_step is None:\n            first_step = 0  # LSODA value for automatic selection.\n        else:\n            first_step = validate_first_step(first_step, t0, t_bound)\n\n        first_step *= self.direction\n\n        if max_step == np.inf:\n            max_step = 0  # LSODA value for infinity.\n        elif max_step <= 0:\n            raise ValueError(\"`max_step` must be positive.\")\n\n        if min_step < 0:\n            raise ValueError(\"`min_step` must be nonnegative.\")\n\n        rtol, atol = validate_tol(rtol, atol, self.n)\n\n        solver = ode(self.fun, jac)\n        solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,\n                              min_step=min_step, first_step=first_step,\n                              lband=lband, uband=uband)\n        solver.set_initial_value(y0, t0)\n\n        # Inject t_bound into rwork array as needed for itask=5.\n        solver._integrator.rwork[0] = self.t_bound\n        solver._integrator.call_args[4] = solver._integrator.rwork\n\n        self._lsoda_solver = solver\n\n    def _step_impl(self):\n        solver = self._lsoda_solver\n        integrator = solver._integrator\n\n        # From lsoda.step and lsoda.integrate itask=5 means take a single\n        # step and do not go past t_bound.\n        itask = integrator.call_args[2]\n        integrator.call_args[2] = 5\n        solver._y, solver.t = integrator.run(\n            solver.f, solver.jac or (lambda: None), solver._y, solver.t,\n            self.t_bound, solver.f_params, solver.jac_params)\n        integrator.call_args[2] = itask\n\n        if solver.successful():\n            self.t = solver.t\n            self.y = solver._y\n            # From LSODA Fortran source njev is equal to nlu.\n            self.njev = integrator.iwork[12]\n            self.nlu = integrator.iwork[12]\n            return True, None\n        else:\n            return False, 'Unexpected istate in LSODA.'\n\n    def _dense_output_impl(self):\n        iwork = self._lsoda_solver._integrator.iwork\n        rwork = self._lsoda_solver._integrator.rwork\n\n        order = iwork[14]\n        h = rwork[11]\n        yh = np.reshape(rwork[20:20 + (order + 1) * self.n],\n                        (self.n, order + 1), order='F').copy()\n\n        return LsodaDenseOutput(self.t_old, self.t, h, order, yh)\n\n\nclass LsodaDenseOutput(DenseOutput):\n    def __init__(self, t_old, t, h, order, yh):\n        super(LsodaDenseOutput, self).__init__(t_old, t)\n        self.h = h\n        self.yh = yh\n        self.p = np.arange(order + 1)\n\n    def _call_impl(self, t):\n        if t.ndim == 0:\n            x = ((t - self.t) / self.h) ** self.p\n        else:\n            x = ((t - self.t) / self.h) ** self.p[:, None]\n\n        return np.dot(self.yh, x)\n", "path": "scipy/integrate/_ivp/lsoda.py"}]}