A Practical Guide to scipy.integrate.solve_ivp

The differential equations in my physics homework all had neat closed-form answers, so the first one from an actual project threw me. A cooling model with awkward constants, no textbook solution, and a deadline.

solve_ivp from scipy.integrate exists for that situation, taking an equation, a time span, and a starting value, and handing back numbers you can plot or feed into the next calculation.

What an initial value problem is

An initial value problem has two parts, a rule for how something changes and a known starting state. Newton’s law of cooling is the classic example, a hot drink losing heat in proportion to how much warmer it is than the room.

Written out, that is dT/dt = -k(T – T_env) with a start like T(0) = 90, and solve_ivp turns the pair into T at any later time.

The function signature that matters

Three arguments are required and their shapes cause most first-run errors. Everything else is optional tuning.

from scipy.integrate import solve_ivp

solution = solve_ivp(
    fun,      # your derivative function, called as fun(t, y)
    t_span,   # tuple: (start_time, end_time)
    y0,       # list or array of starting values, even for one equation
)

The derivative function receives time first and state second. y0 must be a list or array even for a single equation, because solve_ivp treats every problem as a system.

Passing a bare number like 90 instead of [90] raises ValueError with the message “y0 must be 1-dimensional.” That exact error greeted me on my first attempt.

Solving the cooling problem

Here is the full model, a 90 degree coffee in a 22 degree room with cooling constant 0.3, chosen because the exact answer is known and the solver can be checked against it.

import numpy as np
from scipy.integrate import solve_ivp

def cooling(t, T):
    return -0.3 * (T - 22)

solution = solve_ivp(cooling, (0, 10), [90])

print(solution.y[0, -1])
print(22 + (90 - 22) * np.exp(-0.3 * 10))

Output:

25.391747438694296
25.38552674282866

After ten minutes the numerical answer is 25.3917 against an exact 25.3855, a gap of about 0.006 degrees. The solver got there in only six time steps because it stretches the step size wherever the solution is smooth.

Reading the result object

solve_ivp returns a bunch object rather than a plain array, and four fields cover almost everything you will do with it.

The row-per-equation layout means solution.y[0] is the full time series of your first variable, and y[0, -1] is its final value.

Getting values at the times you actually want

The solver’s automatically chosen time points are efficient but irregular, which is awkward for plots and tables. The t_eval argument requests output at specific times without changing how the integration runs.

import numpy as np
from scipy.integrate import solve_ivp

def cooling(t, T):
    return -0.3 * (T - 22)

times = np.linspace(0, 10, 5)
solution = solve_ivp(cooling, (0, 10), [90], t_eval=times)

print(solution.t)
print(np.round(solution.y[0], 2))

Output:

[ 0.   2.5  5.   7.5 10. ]
[90.   54.1  37.19 29.18 25.39]

Five evenly spaced readings tell the whole story, a fast drop early and a slow crawl toward room temperature at the end, which is exponential decay doing what it always does.

Second-order equations become systems

solve_ivp only understands first derivatives, and plenty of interesting equations have second ones. The standard trick rewrites x′′ as a pair of first-order equations by treating position and velocity as separate variables.

A damped oscillator x′′ + 0.5x′ + 4x = 0 becomes two rules, position changing at the rate of velocity and velocity changing under the spring and damping forces.

import numpy as np
from scipy.integrate import solve_ivp

def oscillator(t, state):
    x, v = state
    return [v, -0.5 * v - 4 * x]

solution = solve_ivp(
    oscillator, (0, 10), [1, 0],
    t_eval=np.linspace(0, 10, 6),
)

print(np.round(solution.y[0], 4))

Output:

[ 1.     -0.4664  0.0151  0.1589 -0.1358  0.053 ]

The position swings negative, comes back, and shrinks toward zero as damping drains the energy. y0 now holds two numbers, the starting position 1 and starting velocity 0, matching the two equations.

Passing constants with args

Hardcoding k and the room temperature inside the function gets old around the third time you change them, and the args parameter is the way out, forwarding extra values to every call of your derivative function.

from scipy.integrate import solve_ivp

def cooling(t, T, k, T_env):
    return -k * (T - T_env)

solution = solve_ivp(cooling, (0, 10), [90], args=(0.3, 22))
print(solution.y[0, -1])

Same answer as before, 25.3917, with the constants now living at the call site. Sweeping k across a range becomes a loop over args instead of an edit-and-rerun cycle.

Stopping when something happens with events

For a ball thrown upward, the number I actually want is the landing time. Integrating to some arbitrary end point and then digging through arrays for the zero crossing is backwards.

An event function returns a value that crosses zero at the moment you care about. Mark it terminal and the solver stops right there.

from scipy.integrate import solve_ivp

def ball(t, state):
    y, v = state
    return [v, -9.81]

def hit_ground(t, state):
    return state[0]

hit_ground.terminal = True
hit_ground.direction = -1

solution = solve_ivp(ball, (0, 20), [0, 12], events=hit_ground)
print(solution.t_events[0])

A ball thrown up at 12 m/s lands at t = 2.4465 seconds, and the kinematics formula 2v/g gives 2.4465 as well. The direction attribute set to -1 means only downward crossings count, so the launch moment at height zero is ignored.

When the default solver struggles

The default method, RK45, handles most smooth problems well. Stiff equations make it crawl. Stiffness means some part of the solution decays enormously faster than the rest, and RK45 shrinks its steps to chase the fast part long after it has died away.

import numpy as np
from scipy.integrate import solve_ivp

def stiff(t, y):
    return -1000 * (y - np.cos(t))

fast = solve_ivp(stiff, (0, 2), [0], method="Radau")
slow = solve_ivp(stiff, (0, 2), [0])

print(slow.nfev, "function calls with RK45")
print(fast.nfev, "function calls with Radau")

Output:

4238 function calls with RK45
100 function calls with Radau

Same equation, same accuracy target, and the implicit Radau method needs 42 times fewer function evaluations. When a run feels mysteriously slow, or nfev looks huge, switching method to Radau or BDF is the first thing to try.

Conclusion

solve_ivp turns a differential equation into arrays, with one function call standing between the model and the numbers. The shape rules catch beginners first, y0 as a list and fun taking t before y.

From there, t_eval controls the output grid, args carries constants, events catches moments, and method rescues stiff problems. That handful of arguments has covered every initial value problem I have hit since the cooling model that started this.