How to Imitate Ode45() Function in Python

The ode45() function in MATLAB solves initial value problems for ordinary differential equations (ODEs) with an adaptive Runge-Kutta method. Python has no function literally named ode45(), but SciPy’s solve_ivp() covers the same workflow: you describe the equation, give the solver a time interval and a starting state, and read the solution from the object it returns.

Use SciPy’s solve_ivp() for MATLAB’s ode45()

The solve_ivp() function in scipy.integrate is the closest practical analogue to ode45(). Its default method, RK45, is an explicit Runge-Kutta formula of order 5(4) with automatic step-size control — the same algorithm family that ode45() uses in MATLAB, so the solver takes larger steps where the solution changes slowly and smaller steps where it changes quickly.

It is not a drop-in replacement: the argument names differ and the result comes back as a single object rather than a pair of t and y arrays. Still, the overall pattern — a derivative function, a time span, and an initial state — maps directly between the two. Install SciPy with pip install scipy if it is not already available.

Define the ODE and call solve_ivp()

The essential call is solve_ivp(fun, t_span, y0, method="RK45", t_eval=None, args=None, ...). The main inputs are:

• fun: a callable fun(t, y) that returns the derivative dy/dt for the current time t and state y. For a single equation it returns one value; for a system it returns an array with one derivative per state component.
• t_span: a pair (t_start, t_end) giving the start and end of the integration interval.
• y0: the state at t_start; its length sets the number of equations in the system.
• method (optional): the integration scheme, 'RK45' by default.
• args (optional): extra values passed to the function as fun(t, y, *args).
• rtol and atol (optional): relative and absolute tolerances controlling the adaptive step size.

This example solves dy/dt = 2y starting from y(0) = 1 over the interval [0, 2]. The call returns a result object; printing its t and y attributes shows the time points the solver chose and the state at each of them:

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


result = solve_ivp(growth, (0, 2), [1])
print("t:", result.t)
print("y:", result.y[0])

t: [0.         0.07580098 0.56590702 1.18757739 1.85907272 2.        ]
y: [ 1.          1.16369695  3.10133569 10.75272756 41.18553629 54.59497095]

The spacing of the reported times is uneven because solve_ivp() adapts each step to the local behavior of the solution rather than marching at a fixed increment.

Read the solver result

solve_ivp() returns a scipy.integrate.OdeResult object. The attributes you will use most are:

• t: a NumPy array of the reported time points. Without t_eval, these are the solver’s internal adaptive steps.
• y: a 2-D NumPy array of shape (n, len(t)); row i holds component i of the state at each reported time, so y[:, k] is the state at t[k]. For a single equation, y[0] is the solution curve.
• sol: a continuous OdeSolution interpolant you can call as result.sol(t) to evaluate the solution at arbitrary times inside the interval. It is only populated when you pass dense_output=True; otherwise it is None.
• success and message: a Boolean for whether the integration finished and a human-readable description of how it ended.
• status: 0 when the interval was completed, 1 when a termination event occurred, and -1 when the solver failed to make progress.

Checking success (or status) before using y is a good habit, since a failed run can still return partial arrays.

Choose output times with t_eval

The adaptive step times rarely match the times you want to report. Pass a sorted array of times to t_eval and solve_ivp() returns the solution at exactly those times instead of its internal steps:

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


times = [0.5, 1.0, 1.5, 2.0]
result = solve_ivp(growth, (0, 2), [1], t_eval=times)
print("t:", result.t)
print("y:", result.y[0])

t: [0.5 1.  1.5 2. ]
y: [ 2.71823524  7.38683161 20.07910026 54.59497095]

Now t contains only the requested times and y holds the interpolated state at each one. Every value in t_eval must lie inside t_span; asking for a time outside the integration interval stops the call with an error stating that the values must be within t_span. Endpoints are allowed, so t_eval may include both t_start and t_end.