What’s actually slowing this PC down?
Pick the symptom - the matching free tool is one click away.
To solve a differential equation with scipy.integrate.odeint, write a function that returns the derivatives, pass it an initial state and an array of output times, and read the answer from the returned array. odeint still works, but SciPy’s own reference says: “For new code, use scipy.integrate.solve_ivp to solve a differential equation.” This guide covers both: the odeint pattern you need for existing code, and the solve_ivp equivalent for anything new.
Version context: the odeint reference page checked here is labelled v1.11.4, while the integration tutorial and solve_ivp reference are v1.18.0. The SciPy homepage lists 1.18.1 (released 2026-08-21). Check your installed version with import scipy; print(scipy.__version__).
As an Amazon Associate I earn from qualifying purchases.
The minimal odeint pattern
odeint solves initial-value problems for systems of first-order ODEs using LSODA from the ODEPACK library, which handles both stiff and non-stiff problems. Here is the decay equation dy/dt = −k·y with y(0) = 1:
import numpy as np
from scipy.integrate import odeint
def decay(y, t, k): # default order: (y, t, ...)
return -k * y
t = np.linspace(0.0, 5.0, 101)
y0 = 1.0
k = 0.7
solution = odeint(decay, y0, t, args=(k,))
y = solution[:, 0] # 1-D vector for plottingThis example is built from the documented signature rather than taken from SciPy’s docs. Its result should approximate y(t) = e−0.7t, which makes it a handy sanity check.
#1 Best Overall
What each argument does
- func: returns dy/dt. By default it is called as
func(y, t, ...). Passtfirst=Trueto usefunc(t, y, ...)instead. - y0: the initial state, a scalar or a sequence with one value per state variable.
- t: the times at which you want output. The sequence must be monotonically increasing or decreasing; repeated values are allowed. The first entry is the initial time.
- args: a tuple of extra parameters forwarded to
func. A single parameter still needs the trailing comma:args=(k,).
Reading the result
The output has shape (len(t), len(y0)): one row per requested time, one column per state variable, with the initial state in row zero. For a scalar problem, take solution[:, 0]. For a system, solution[:, i] is the trajectory of variable i.
Turning higher-order equations into first-order systems
Solvers accept only first-order systems. For a second-order equation x″ = g(x, x′, t), introduce y[0] = x and y[1] = x′. The derivative function returns [y[1], g(y[0], y[1], t)], and y0 holds both x(0) and x′(0). SciPy’s tutorial uses the same conversion principle for its second-order example. Higher orders work the same way, adding one state per extra derivative.
A damped oscillator, x″ + c·x′ + w²·x = 0:
def oscillator(y, t, c, w):
x, v = y
return [v, -c * v - w**2 * x]
t = np.linspace(0, 20, 400)
sol = odeint(oscillator, [1.0, 0.0], t, args=(0.3, 2.0))
x = sol[:, 0] # position
v = sol[:, 1] # velocityThe same problem with solve_ivp
solve_ivp is SciPy’s recommended interface. The differences are small but easy to trip over:
Recommended Free Tools
from scipy.integrate import solve_ivp
def decay(t, y, k): # note: (t, y)
return -k * y
res = solve_ivp(decay, (0.0, 5.0), [1.0], args=(0.7,),
t_eval=np.linspace(0.0, 5.0, 101))
t_out = res.t
y_out = res.y[0] # states on rows, times on columns- Pass the interval
(t0, tf)ast_span; the solver chooses its own steps. Uset_evalonly if you want output at specific times. - The initial state must be array-like, so use
[1.0]even for one variable. - The result is an object. Check
res.successandres.messagebefore trustingres.y.
odeint versus solve_ivp
| Decision axis | odeint | solve_ivp |
|---|---|---|
| SciPy guidance | Reasonable for existing code; the reference recommends solve_ivp for new code |
Recommended for new code |
| Callback order | func(y, t, ...) by default; tfirst=True switches |
fun(t, y) |
| Time input | Array of requested output times | Interval t_span, optional t_eval |
| Result layout | Array shaped (len(t), len(y0)) |
Result object; y has states on rows, times on columns |
| Solvers | LSODA | RK45 (default), RK23, DOP853, Radau, BDF, LSODA |
| Extras highlighted in docs | Optional Jacobian, diagnostics, banded-Jacobian controls | Events, dense output, solver choice, status information |
Moving old code across usually means swapping the argument order, replacing the time array with t_span plus t_eval, and transposing the output (res.y.T matches the odeint layout).
Choosing a solver and handling stiffness
SciPy recommends explicit Runge–Kutta methods (RK45, RK23, DOP853) for non-stiff problems and implicit methods (Radau, BDF) for stiff ones. For solve_ivp the documentation advises: “If not sure, first try to run ‘RK45’.” If it needs an unusually large number of iterations or fails to integrate, the problem is probably stiff, so switch to Radau or BDF. LSODA is also available there and is what odeint uses.
With odeint, you can supply a Jacobian through Dfun when you have one. If the Jacobian is banded, ml and mu give its lower and upper bandwidths. SciPy’s tutorial demonstrates the payoff on a 5,000-state Gray–Scott example: 25.2 seconds per loop without band information versus 191 milliseconds per loop with ml=2 and mu=2 (tutorial v1.18.0). That is one example’s timing, not a general speed-up figure.
Tolerances: what they do and do not guarantee
Both interfaces expose rtol and atol, which control the local error estimate at each step. Set them with the scale of your variables in mind: an absolute tolerance far larger than a state’s typical magnitude lets that variable drift unchecked. solve_ivp accepts an array of per-component atol values for this reason.
Local error control does not prove global accuracy. Validate results against a known analytical solution, conserved quantities, or a rerun with tighter tolerances. SciPy’s tutorial does this with the Airy function, showing better agreement after tightening tolerances.
Common mistakes
- Swapped arguments. A function written as
f(t, y)receives wrong values under defaultodeint. Rewrite it asf(y, t)or passtfirst=True. Often nothing errors; the numbers are simply wrong. - Passing an interval to odeint.
odeintwants the output times, not(t0, tf). - Wrong indexing.
odeintoutput is time-by-state;solve_ivp‘sres.yis state-by-time. - Forgetting to reduce the order. Add a state for each derivative below the highest.
- Returning the wrong shape. The derivative function must return one value per state variable, in the same order as
y0.
The Bottom Line
Use solve_ivp for new projects and keep odeint for existing code, where it remains documented and functional. In both cases, reduce to first order, match the argument order to the solver, and check the result against something you know.
Quick Recap
Product prices and availability are accurate as of the date/time indicated and are subject to change. Any price and availability information displayed on Amazon at the time of purchase will apply.
