MATLAB Function Reference 
Solve initial value problems for ordinary differential equations (ODEs)
Syntax
[T,Y] = solver(odefun,tspan,y0) [T,Y] = solver(odefun,tspan,y0,options) [T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options) sol = solver(odefun,[t0 tf],y0...)
where solver
is one of ode45
, ode23
, ode113
, ode15s
, ode23s
, ode23t
, or ode23tb
.
Arguments
The following table describes the input arguments to the solvers.

A function handle that evaluates the right side of the differential equations. See Function Handles in the MATLAB Programming documentation for more information. All solvers solve systems of equations in the form or problems that involve a mass matrix, . The ode23s solver can solve only equations with constant mass matrices. ode15s and ode23t can solve problems with a mass matrix that is singular, i.e., differentialalgebraic equations (DAEs). 
tspan 
A vector specifying the interval of integration, [t0,tf] . The solver imposes the initial conditions at tspan(1) , and integrates from tspan(1) to tspan(end) . To obtain solutions at specific times (all increasing or all decreasing), use tspan = [t0,t1,...,tf] .For tspan vectors with two elements [t0 tf] , the solver returns the solution evaluated at every integration step. For tspan vectors with more than two elements, the solver returns solutions evaluated at the given time points. The time values must be in order, either increasing or decreasing. 

Specifying tspan with more than two elements does not affect the internal time steps that the solver uses to traverse the interval from tspan(1) to tspan(end) . All solvers in the ODE suite obtain output values by means of continuous extensions of the basic formulas. Although a solver does not necessarily step precisely to a time point specified in tspan , the solutions produced at the specified time points are of the same order of accuracy as the solutions computed at the internal time points. Specifying tspan with more than two elements has little effect on the efficiency of computation, but for large systems, affects memory management. 
y0 
A vector of initial conditions. 

Structure of optional parameters that change the default integration properties. This is the fourth input argument. You can create options using the odeset function. See odeset for details. 
The following table lists the output arguments for the solvers.
T 
Column vector of time points 
Y 
Solution array. Each row in y corresponds to the solution at a time returned in the corresponding row of t . 
Description
[T,Y] =
with solver
(odefun,tspan,y0)
tspan = [t0 tf]
integrates the system of differential equations from time t0
to tf
with initial conditions y0
. odefun
is a function handle. See Function Handles in the MATLAB Programming documentation for more information. Function f = odefun(t,y)
, for a scalar t
and a column vector y
, must return a column vector f
corresponding to . Each row in the solution array Y
corresponds to a time returned in column vector T
. To obtain solutions at the specific times t0
, t1,...,tf
(all increasing or all decreasing), use tspan = [t0,t1,...,tf]
.
Parameterizing Functions Called by Function Functions, in the MATLAB mathematics documentation, 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 in solver
(odefun,tspan,y0,options)
options
, an argument created with the odeset
function. Commonly used properties include a scalar relative error tolerance RelTol
(1e3
by default) and a vector of absolute error tolerances AbsTol
(all component
s are 1e6
by default). See odeset
for details.
[T,Y,TE,YE,IE] =
solves as above while also finding where functions of , 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 solver
(odefun,tspan,y0,options)
'Events'
property to a function, e.g., events
or @events
, and
creating a function [value
,isterminal
,direction
] = events
(t
,y
). For the i
th event function in events
:
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 0
otherwise.
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 TE
, YE
, and IE
return, respectively, the time at which an event occurs, the solution at the time of the event, and the index i
of the event function that vanishes.
sol =
returns a structure that you can use with solver
(odefun,[t0 tf],y0...)
deval
to evaluate the solution at any point on the interval [t0,tf]
. You must pass odefun
as a function handle. The structure sol
always includes these fields:
sol.x 
Steps chosen by the solver. 
sol.y 
Each column sol.y(:,i) contains the solution at sol.x(i) . 
sol.solver 
Solver name. 
If you specify the Events
option and events are detected, sol
also includes these fields:
If you specify an output function as the value of the OutputFcn
property, the solver calls it with the computed solution after each time step. Four output functions are provided: odeplot
, odephas2
, odephas3
, odeprint
. When you call the solver with no output arguments, it calls the default odeplot
to plot the solution as it is computed. odephas2
and odephas3
produce two and threedimensional phase plane plots, respectively. odeprint
displays 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 of OutputSel
to [1,3]
, the solver plots solution components 1 and 3 as they are computed.
For the stiff solvers ode15s
, ode23s
, ode23t
, and ode23tb
, the Jacobian matrix is critical to reliability and efficiency. Use odeset
to set Jacobian
to @FJAC
if FJAC(T,Y)
returns the Jacobian or to the matrix if the Jacobian is constant. If the Jacobian
property is not set (the default), is approximated by finite differences. Set the Vectorized
property 'on
' if the ODE function is coded so that odefun
(T
,[Y1,Y2 ...
]) returns [odefun
(T
,Y1
),odefun
(T
,Y2
) ...
]. If is a sparse matrix, set the JPattern
property to the sparsity pattern of , i.e., a sparse matrix S
with S(i,j)
= 1 if the i
th component of depends on the j
th component of , and 0 otherwise.
The solvers of the ODE suite can solve problems of the form , with time and statedependent mass matrix . (The ode23s
solver 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 odeset
to set the Mass
property to @MASS
. If the mass matrix is constant, the matrix should be used as the value of the Mass
property. Problems with statedependent mass matrices are more difficult:
MASS
is to be called with one input argument, t
, set the MStateDependence
property to 'none
'.
MStateDependence
to 'weak
' (the default) and otherwise, to 'strong
'. In either case, the function MASS
is called with the two arguments (t
,y
).
If there are many differential equations, it is important to exploit sparsity:
JPattern
property or a sparse using the Jacobian
property.
MvPattern
to a sparse matrix S
with S(i,j) = 1
if for any k
, the (i,k
) component of depends on component j
of , and 0
otherwise.
If the mass matrix is singular, then is a differential algebraic equation. DAEs have solutions only when is consistent, that is, if there is a vector such that . The ode15s
and ode23t
solvers 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 MassSingular
property to 'yes'
, 'no'
, or 'maybe'
. The default value of 'maybe'
causes the solver to test whether the problem is a DAE. You can provide yp0
as the value of the InitialSlope
property. The default is the zero vector. If a problem is a DAE, and y0
and yp0
are 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 is a diagonal matrix (a semiexplicit DAE).
The algorithms used in the ODE solvers vary according to order of accuracy [6] and the type of systems (stiff or nonstiff) they are designed to solve. See Algorithms for more details.
Options
Different solvers accept different parameters in the options list. For more information, see odeset
and Changing ODE Integration Properties in the MATLAB documentation.
Examples
Example 1. 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 rigid
containing the equations
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 [0 12]
with an initial condition vector [0 1 1]
at time 0
.
options = odeset('RelTol',1e4,'AbsTol',[1e4 1e4 1e5]); [T,Y] = ode45(@rigid,[0 12],[0 1 1],options);
Plotting the columns of the returned array Y
versus T
shows the solution
Example 2. 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 vdp1000
containing the equations
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 tolerances (1e3
and 1e6
, respectively) and solve on a time interval of [0 3000]
with initial condition vector [2 0]
at time 0
.
Plotting the first column of the returned matrix Y
versus T
shows the solution
Algorithms
ode45
is based on an explicit RungeKutta (4,5) formula, the DormandPrince pair. It is a onestep solver  in computing y(t
_{n})
, it needs only the solution at the immediately preceding time point, y(t
_{n1})
. In general, ode45
is the best function to apply as a "first try" for most problems. [3]
ode23
is an implementation of an explicit RungeKutta (2,3) pair of Bogacki and Shampine. It may be more efficient than ode45
at crude tolerances and in the presence of moderate stiffness. Like ode45
, ode23
is a onestep solver. [2]
ode113
is a variable order AdamsBashforthMoulton PECE solver. It may be more efficient than ode45
at 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. [7]
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 ode113
, ode15s
is a multistep solver. Try ode15s
when ode45
fails, or is very inefficient, and you suspect that the problem is stiff, or when solving a differentialalgebraic problem. [9], [10]
ode23s
is based on a modified Rosenbrock formula of order 2. Because it is a onestep solver, it may be more efficient than ode15s
at crude tolerances. It can solve some kinds of stiff problems for which ode15s
is not effective. [9]
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 numerical damping. ode23t
can solve DAEs. [10]
ode23tb
is an implementation of TRBDF2, an implicit RungeKutta 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. [8], [1]
See Also
deval
, ode15i
, odeget
, odeset
, function_handle
(@
)
References
[1] 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 436451.
[2] Bogacki, P. and L. F. Shampine, "A 3(2) pair of RungeKutta formulas," Appl. Math. Letters, Vol. 2, 1989, pp 19.
[3] Dormand, J. R. and P. J. Prince, "A family of embedded RungeKutta formulae," J. Comp. Appl. Math., Vol. 6, 1980, pp 1926.
[4] Forsythe, G. , M. Malcolm, and C. Moler, Computer Methods for Mathematical Computations, PrenticeHall, New Jersey, 1977.
[5] Kahaner, D. , C. Moler, and S. Nash, Numerical Methods and Software, PrenticeHall, New Jersey, 1989.
[6] Shampine, L. F. , Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New York, 1994.
[7] Shampine, L. F. and M. K. Gordon, Computer Solution of Ordinary Differential Equations: the Initial Value Problem, W. H. Freeman, San Francisco, 1975.
[8] Shampine, L. F. and M. E. Hosea, "Analysis and Implementation of TRBDF2," Applied Numerical Mathematics 20, 1996.
[9] Shampine, L. F. and M. W. Reichelt, "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp 122.
[10] Shampine, L. F., M. W. Reichelt, and J.A. Kierzenka, "Solving Index1 DAEs in MATLAB and Simulink," SIAM Review, Vol. 41, 1999, pp 538552.
ode15i  odefile 
© 19942005 The MathWorks, Inc.