Solve nonstiff differential equations — medium order method
tspan = [t0 tf], integrates the system of
differential equations from
y0. Each row in the solution
y corresponds to a value returned in column
All MATLAB® ODE solvers can solve systems of equations of
the form ,
or problems that involve a mass matrix, .
The solvers all use similar syntaxes. The
only can solve problems with a mass matrix if the mass matrix is constant.
solve problems with a mass matrix that is singular, known as differential-algebraic
equations (DAEs). Specify the mass matrix using the
ode45 is a versatile ODE solver and is the
first solver you should try for most problems. However, if the problem
is stiff or requires high accuracy, then there are other ODE solvers
that might be better suited to the problem. See Choose an ODE Solver for
finds where functions of (t,y),
called event functions, are zero. In the output,
the time of the event,
ye is the solution at the
time of the event, and
ie is the index of the triggered
For each event function, 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, such as
and creating a corresponding function: [
For more information, see ODE Event Location.
Simple ODEs that have a single solution component can be specified as an anonymous function in the call to the solver. The anonymous function must accept two inputs
(t,y) even if one of the inputs is not used.
Solve the ODE
Use a time interval of
[0,5] and the initial condition
y0 = 0.
tspan = [0 5]; y0 = 0; [t,y] = ode45(@(t,y) 2*t, tspan, y0);
Plot the solution.
The van der Pol equation is a second order ODE
where is a scalar parameter. Rewrite this equation as a system of first-order ODEs by making the substitution . The resulting system of first-order ODEs is
The function file
vdp1.m represents the van der Pol equation using . The variables and are the entries
y(2) of a two-element vector,
function dydt = vdp1(t,y) %VDP1 Evaluate the van der Pol ODEs for mu = 1 % % See also ODE113, ODE23, ODE45. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
Solve the ODE using the
ode45 function on the time interval
[0 20] with initial values
[2 0]. The resulting output is a column vector of time points
t and a solution array
y. Each row in
y corresponds to a time returned in the corresponding row of
t. The first column of
y corresponds to , and the second column to .
[t,y] = ode45(@vdp1,[0 20],[2; 0]);
Plot the solutions for and against
plot(t,y(:,1),'-o',t,y(:,2),'-o') title('Solution of van der Pol Equation (\mu = 1) with ODE45'); xlabel('Time t'); ylabel('Solution y'); legend('y_1','y_2')
ode45 only works with functions that use two input arguments,
y. However, you can pass in extra parameters by defining them outside the function and passing them in when you specify the function handle.
Solve the ODE
Rewriting the equation as a first-order system yields
odefcn.m represents this system of equations as a function that accepts four input arguments:
function dydt = odefcn(t,y,A,B) dydt = zeros(2,1); dydt(1) = y(2); dydt(2) = (A/B)*t.*y(1);
Solve the ODE using
ode45. Specify the function handle such that it passes in the predefined values for
A = 1; B = 2; tspan = [0 5]; y0 = [0 0.01]; [t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);
Plot the results.
For simple ODE systems with one equation, you can specify
y0 as a vector containing multiple initial conditions. This technique creates a system of independent equations through scalar expansion, one for each initial value, and
ode45 solves the system to produce results for each initial value.
Create an anonymous function to represent the equation . The function must accept two inputs for
yprime = @(t,y) -2*y + 2*cos(t).*sin(2*t);
Create a vector of different initial conditions in the range .
y0 = -5:5;
Solve the equation for each initial condition over the time interval using
tspan = [0 3]; [t,y] = ode45(yprime,tspan,y0);
Plot the results.
plot(t,y) grid on xlabel('t') ylabel('y') title('Solutions of y'' = -2y + 2 cos(t) sin(2t), y(0) = -5,-4,...,4,5','interpreter','latex')
This technique is useful for solving simple ODEs with several initial conditions. However, the technique also has some tradeoffs:
You cannot solve systems of equations with multiple initial conditions. The technique only works when solving one equation with multiple initial conditions.
The time step chosen by the solver at each step is based on the equation in the system that needs to take the smallest step. This means the solver can take small steps to satisfy the equation for one initial condition, but the other equations, if solved on their own, would use different step sizes. Despite this, solving for multiple initial conditions at the same time is generally faster than solving the equations separately using a
Consider the following ODE with time-dependent parameters
The initial condition is . The function
f(t) is defined by the n-by-1 vector
f evaluated at times
ft. The function
g(t) is defined by the m-by-1 vector
g evaluated at times
Create the vectors
ft = linspace(0,5,25); f = ft.^2 - ft - 3; gt = linspace(1,6,25); g = 3*sin(gt-0.25);
Write a function named
myode that interpolates
g to obtain the value of the time-dependent terms at the specified time. Save the function in your current folder to run the rest of the example.
myode function accepts extra input arguments to evaluate the ODE at each time step, but
ode45 only uses the first two input arguments
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
Solve the equation over the time interval
[1 5] using
ode45. Specify the function using a function handle so that
ode45 uses only the first two input arguments of
myode. Also, loosen the error thresholds using
tspan = [1 5]; ic = 1; opts = odeset('RelTol',1e-2,'AbsTol',1e-4); [t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);
Plot the solution,
y, as a function of the time points,
The van der Pol equation is a second order ODE
Solve the van der Pol equation with using
ode45. The function
vdp1.m ships with MATLAB® and encodes the equations. Specify a single output to return a structure containing information about the solution, such as the solver and evaluation points.
tspan = [0 20]; y0 = [2 0]; sol = ode45(@vdp1,tspan,y0)
sol = struct with fields: solver: 'ode45' extdata: [1x1 struct] x: [1x60 double] y: [2x60 double] stats: [1x1 struct] idata: [1x1 struct]
linspace to generate 250 points in the interval
[0 20]. Evaluate the solution at these points using
x = linspace(0,20,250); y = deval(sol,x);
Plot the first component of the solution.
Extend the solution to using
odextend and add the result to the original plot.
sol_new = odextend(sol,@vdp1,35); x = linspace(20,35,350); y = deval(sol_new,x); hold on plot(x,y(1,:),'r')
ode45 is based on an explicit Runge-Kutta
(4,5) formula, the Dormand-Prince pair. It is a single-step solver
– in computing
it needs only the solution at the immediately preceding time point,
y(tn-1) , .
 Dormand, J. R. and P. J. Prince, “A family of embedded Runge-Kutta formulae,” J. Comp. Appl. Math., Vol. 6, 1980, pp. 19–26.
 Shampine, L. F. and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.
Usage notes and limitations:
arguments must be constant.
Code generation does not support a constant mass matrix in the options structure. Provide a mass matrix as a function.
You must provide at least the two output arguments
Input types must be homogeneous—all double or all single.
Variable-sizing support must be enabled. Code generation
requires dynamic memory allocation when
two elements or you use event functions.