Solving differential equations with solve_ivp

Integrate an initial value problem, pick a solver that matches the physics, detect events, and recognise the cost of an unnecessarily explicit method.

The shape of the problem

import numpy as np
from scipy.integrate import solve_ivp

def decay(t, y):
    """dN/dt = -k N, written for a state vector."""
    k = 0.4
    return [-k * y[0]]

sol = solve_ivp(decay, t_span=(0, 10), y0=[100.0], method="RK45",
                t_eval=np.linspace(0, 10, 200), rtol=1e-8, atol=1e-10)

print(sol.t.shape, sol.y.shape)     # (200,) and (1, 200)
print(sol.success, sol.nfev)

The function signature is f(t, y) and must return the derivative with the same shape as y, even for a scalar problem.
y0 is a sequence; sol.y has shape (n_states, n_times), not the other way round.
Without t_eval the solver returns whichever steps it took; with it you get values at exactly those points, interpolated from the internal solution.

Stiff versus non-stiff

Method	Type	Use for
`RK45`	Explicit, default	Smooth, non-stiff problems; cheap per step
`DOP853`	Explicit, high order	Very smooth problems needing high accuracy
`RK23`	Explicit, low order	Tolerant accuracy, quick exploration
`Radau`	Implicit	Stiff systems needing tight tolerances
`BDF`	Implicit, multistep	Stiff systems with large, sparse Jacobians
`LSODA`	Switches automatically	When you do not know whether it is stiff

def robertson(t, y):
    k1, k2, k3 = 0.04, 3e7, 1e4
    return [-k1 * y[0] + k3 * y[1] * y[2],
             k1 * y[0] - k2 * y[1] ** 2 - k3 * y[1] * y[2],
             k2 * y[1] ** 2]

fast = solve_ivp(robertson, (0, 1e5), [1, 0, 0], method="RK45")
slow = solve_ivp(robertson, (0, 1e5), [1, 0, 0], method="BDF")
print(fast.nfev, slow.nfev)     # the explicit solver takes orders of magnitude more function calls

Stiffness is not a property of the equation alone but of the ratio of timescales present. If an explicit method forces tiny steps for stability while accuracy alone would allow large ones, the problem is stiff and an implicit method will win.

Events and dense output

def falling(t, y):
    return [y[1], -9.81]

def hits_ground(t, y):
    return y[0]                    # zero when height reaches zero

hits_ground.terminal = True        # stop integrating
hits_ground.direction = -1         # only detect a downward crossing

sol = solve_ivp(falling, (0, 10), [10.0, 0.0], events=hits_ground, dense_output=True)
print(sol.t_events[0])                       # time of impact
print(sol.sol(1.234))                        # value anywhere in the range

from scipy.integrate import solve_bvp
def ode(x, y): return np.vstack([y[1], -y[0]])
def bc(ya, yb): return np.array([ya[0], yb[0] - 1])
x = np.linspace(0, np.pi, 20)
bvp = solve_bvp(ode, bc, x, np.zeros((2, x.size)))
print(bvp.success)

💡

An event function must be continuous and change sign at the event; the solver locates the crossing by root finding over the previous step. A discontinuous event function makes detection unreliable, and a missed event usually means the step size stepped straight over it.

FAQ

The solver returns success=False. What now?

Read sol.message. Usually the maximum number of steps was exceeded because the problem is stiff (switch to BDF or Radau), the tolerances are unrealistically tight, or the derivative returns NaN. An rtol below about 1e-13 is rarely achievable in double precision.

How do I pass extra parameters to the derivative?

Use args=(k1, k2) in solve_ivp, or define the derivative as a closure. Avoid reading module-level globals inside it: the values then cannot change between runs and the function becomes harder to test.

Interpolation, smoothing and integration Linear algebra with scipy.linalg

Last refreshed 2026-09-18.

Solving differential equations with solve_ivp

The shape of the problem

Stiff versus non-stiff

Events and dense output

FAQ

Related