Solve stiff differential equations and DAEs; variable order method
This page contains an overview of the solver functions:
ode23tb. You can call any of these solvers
by substituting the placeholder,
with any of the function names.
The following table describes the input arguments to the solvers.
A function handle that evaluates the right side of the
differential equations. All solvers solve systems of equations in
the form y′ = f(t,y)
or problems that involve a mass matrix, M(t,y)y′ = f(t,y).
A vector specifying the interval of integration,
A vector of initial conditions.
Structure of optional parameters that change the default integration properties. This is the fourth input argument.
[t,y] = solver(odefun,tspan,y0,options)
You can create options using the
The following table lists the output arguments for the solvers.
Column vector of time points.
Solution array. Each row in
The time at which an event occurs.
The solution at the time of the event.
Structure to evaluate the solution.
tspan = [t0 tf] integrates the system of differential
equations y′ = f(t,y)
tf with initial
y0. The first input argument,
is a function handle. The function,
f = odefun(t,y),
for a scalar
t and a column vector
must return a column vector
f corresponding to f(t,y).
Each row in the solution array
Y corresponds to
a time returned in column vector
T. To obtain solutions
at the specific times
increasing or all decreasing), use
tspan = [t0,t1,...,tf].
Parameterizing Functions explains how to provide additional
parameters to the function
fun, if necessary.
[T,Y] = solves
as above with default integration parameters replaced by property values specified
options, an argument created with the
Commonly used properties include a scalar relative error tolerance
default) and a vector of absolute error tolerances
1e-6 by default). If certain components
of the solution must be nonnegative, use the
to set the
NonNegative property to the indices
of these components. See
[T,Y,TE,YE,IE] = solves
as above while also finding where functions of (t,y),
called event functions, are zero. For each event function, you specify
whether the integration is to terminate at a zero and whether the
direction of the zero crossing matters. Do this by setting the
to a function, e.g.,
a function [
ith event function in
value(i) is the value of the function.
isterminal(i) = 1, if the integration
is to terminate at a zero of this event function and
direction(i) = 0 if all zeros are
to be computed (the default),
+1 if only the zeros
where the event function increases, and
-1 if only
the zeros where the event function decreases.
Corresponding entries in
IE return, respectively, the time at which
an event occurs, the solution at the time of the event, and the index
the event function that vanishes.
sol = returns a structure that you can use with
evaluate the solution at any point on the interval
You must pass
odefun as a function handle. The
sol always includes these fields:
Steps chosen by the solver.
If you specify the
Events option and events
sol also includes these fields:
Points at which events, if any, occurred.
Solutions that correspond to events in
Indices into the vector returned by the function specified
If you specify an output function as the value of the
the solver calls it with the computed solution after each time step.
Four output functions are provided:
When you call the solver with no output arguments, it calls the default
plot the solution as it is computed.
two- and three-dimensional phase plane plots, respectively.
the solution components on the screen. By default, the ODE solver
passes all components of the solution to the output function. You
can pass only specific components by providing a vector of indices
as the value of the
OutputSel property. For example,
if you call the solver with no output arguments and set the value
[1,3], the solver
plots solution components 1 and 3 as they are computed.
For the stiff solvers
ode23tb, the Jacobian matrix ∂f/∂y is
critical to reliability and efficiency. Use
the Jacobian ∂f/∂y or
to the matrix ∂f/∂y if
the Jacobian is constant. If the
is not set (the default), ∂f/∂y is
approximated by finite differences. Set the
on' if the ODE function is coded so that
...]) returns [
If ∂f/∂y is a
sparse matrix, set the
JPattern property to the
sparsity pattern of ∂f/∂y,
i.e., a sparse matrix
1 if the
ith component of f(t,y)
depends on the
jth component of y,
and 0 otherwise.
The solvers of the ODE suite can solve problems of the form M(t,y)y′ = f(t,y),
with time- and state-dependent mass matrix M. (The
can solve only equations with constant mass matrices.) If a problem
has a mass matrix, create a function
M = MASS(t,y) that
returns the value of the mass matrix, and use
Mass property to
If the mass matrix is constant, the matrix should be used as the value
Mass property. Problems with state-dependent
mass matrices are more difficult:
If the mass matrix does not depend on the state variable y and
MASS is to be called with one input
t, set the
If the mass matrix depends weakly on y,
MStateDependence to '
(the default); otherwise, set it to '
either case, the function
MASS is called with
the two arguments (
If there are many differential equations, it is important to exploit sparsity:
Return a sparse M(t,y).
Supply the sparsity pattern of ∂f/∂y using
JPattern property or a sparse ∂f/∂y using
For strongly state-dependent M(t,y),
MvPattern to a sparse matrix
= 1 if for any
k, the (
component of M(t,y)
depends on component
j of y,
If the mass matrix M is singular, then M(t,y)y′ = f(t,y)
is a system of differential algebraic equations. DAEs have solutions
only when y0 is consistent,
that is, if there is a vector yp0 such
that M(t0,y0)yp0 = f(t0,y0).
can solve DAEs of index 1 provided that
y0 is sufficiently
close to being consistent. If there is a mass matrix, you can use
odeset to set the
The default value of
'maybe' causes the solver
to test whether the problem is a DAE. You can provide
the value of the
InitialSlope property. The default
is the zero vector. If a problem is a DAE, and
not consistent, the solver treats them as guesses, attempts to compute
consistent values that are close to the guesses, and continues to
solve the problem. When solving DAEs, it is very advantageous to formulate
the problem so that M is a diagonal matrix (a semi-explicit
Order of Accuracy
When to Use
Most of the time. This should be the first solver you try.
For problems with crude error tolerances or for solving moderately stiff problems.
Low to high
For problems with stringent error tolerances or for solving computationally intensive problems.
Low to medium
If using crude error tolerances to solve stiff systems and the mass matrix is constant.
For moderately stiff problems if you need a solution without numerical damping.
If using crude error tolerances to solve stiff systems.
You can use the
An example of a nonstiff system is the system of equations describing the motion of a rigid body without external forces.
To simulate this system, create a function
function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2);
In this example we change the error tolerances using the
odeset command and solve on a time interval
12] with an initial condition vector
[0 1 1] at
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); [T,Y] = ode45(@rigid,[0 12],[0 1 1],options);
Plotting the columns of the returned array
An example of a stiff system is provided by the van der Pol equations in relaxation oscillation. The limit cycle has portions where the solution components change slowly and the problem is quite stiff, alternating with regions of very sharp change where it is not stiff.
To simulate this system, create a function
function dy = vdp1000(t,y) dy = zeros(2,1); % a column vector dy(1) = y(2); dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);
For this problem, we will use the default relative and absolute
and solve on a time interval of
[0 3000] with initial
[2 0] at time
[T,Y] = ode15s(@vdp1000,[0 3000],[2 0]);
Plotting the first column of the returned matrix
This example solves an ordinary differential equation with time-dependent terms.
Consider the following ODE, with time-dependent parameters defined only through the set of data points given in two vectors:
y'(t) + f(t)y(t) = g(t)
y(1) = 1, where the function
f(t)is defined by the
fevaluated at times
ft, and the function
g(t)is defined by the
gevaluated at times
First, define the time-dependent parameters
ft = linspace(0,5,25); % Generate t for f f = ft.^2 - ft - 3; % Generate f(t) gt = linspace(1,6,25); % Generate t for g g = 3*sin(gt-0.25); % Generate g(t)
function dydt = myode(t,y,ft,f,gt,g) f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t dydt = -f.*y + g; % Evaluate ODE at time t
myode.mwithin the MATLAB
ode45function specifying time as the first input argument :
Tspan = [1 5]; % Solve from t=1 to t=5 IC = 1; % y(t=1) = 1 [T Y] = ode45(@(t,y) myode(t,y,ft,f,gt,g),Tspan,IC); % Solve ODE
y(t)as a function of time:
plot(T, Y); title('Plot of y as a function of time'); xlabel('Time'); ylabel('Y(t)');
ode45 is based on an explicit Runge-Kutta
(4,5) formula, the Dormand-Prince pair. It is a one-step solver
– in computing
it needs only the solution at the immediately preceding time point,
y(tn-1). In general,
the best function to apply as a first try for
most problems. 
ode23 is an implementation of an explicit
Runge-Kutta (2,3) pair of Bogacki and Shampine. It may be more efficient
ode45 at crude tolerances and in the presence
of moderate stiffness. Like
a one-step solver. 
ode113 is a variable order Adams-Bashforth-Moulton
PECE solver. It may be more efficient than
stringent tolerances and when the ODE file function is particularly
expensive to evaluate.
ode113 is a multistep solver
— it normally needs the solutions at several preceding time
points to compute the current solution. 
The above algorithms are intended to solve nonstiff systems. If they appear to be unduly slow, try using one of the stiff solvers below.
ode15s is a variable order solver based
on the numerical differentiation formulas (NDFs). Optionally, it uses
the backward differentiation formulas (BDFs, also known as Gear's
method) that are usually less efficient. Like
a multistep solver. Try
or is very inefficient, and you suspect that the problem is stiff,
or when solving a differential-algebraic problem. , 
ode23s is based on a modified Rosenbrock
formula of order 2. Because it is a one-step solver, it may be more
ode15s at crude tolerances. It
can solve some kinds of stiff problems for which
not effective. 
ode23t is an implementation of the trapezoidal
rule using a "free" interpolant. Use this solver if
the problem is only moderately stiff and you need a solution without
ode23t can solve DAEs. 
ode23tb is an implementation of TR-BDF2,
an implicit Runge-Kutta formula with a first stage that is a trapezoidal
rule step and a second stage that is a backward differentiation formula
of order two. By construction, the same iteration matrix is used in
evaluating both stages. Like
ode23s, this solver
may be more efficient than
ode15s at crude tolerances. , 
 Bank, R. E., W. C. Coughran, Jr., W. Fichtner, E. Grosse, D. Rose, and R. Smith, "Transient Simulation of Silicon Devices and Circuits," IEEE Trans. CAD, 4 (1985), pp. 436–451.
 Bogacki, P. and L. F. Shampine, "A 3(2) pair of Runge-Kutta formulas," Appl. Math. Letters, Vol. 2, 1989, pp. 321–325.
 Dormand, J. R. and P. J. Prince, "A family of embedded Runge-Kutta formulae," J. Comp. Appl. Math., Vol. 6, 1980, pp. 19–26.
 Forsythe, G. , M. Malcolm, and C. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, New Jersey, 1977.
 Kahaner, D. , C. Moler, and S. Nash, Numerical Methods and Software, Prentice-Hall, New Jersey, 1989.
 Shampine, L. F. , Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New York, 1994.
 Shampine, L. F. and M. K. Gordon, Computer Solution of Ordinary Differential Equations: the Initial Value Problem, W. H. Freeman, San Francisco, 1975.
 Shampine, L. F. and M. E. Hosea, "Analysis and Implementation of TR-BDF2," Applied Numerical Mathematics 20, 1996.
 Shampine, L. F. and M. W. Reichelt, "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.
 Shampine, L. F., M. W. Reichelt, and J.A. Kierzenka, "Solving Index-1 DAEs in MATLAB and Simulink," SIAM Review, Vol. 41, 1999, pp. 538–552.