This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the page.


Solve fully implicit differential equations — variable order method


  • [t,y] = ode15i(odefun,tspan,y0,yp0)
  • [t,y] = ode15i(odefun,tspan,y0,yp0,options)
  • [t,y,te,ye,ie] = ode15i(odefun,tspan,y0,yp0,options)
  • sol = ode15i(___)



[t,y] = ode15i(odefun,tspan,y0,yp0), where tspan = [t0 tf], integrates the system of differential equations f(t,y,y')=0 from t0 to tf with initial conditions y0 and yp0. Each row in the solution array y corresponds to a value returned in column vector t.


[t,y] = ode15i(odefun,tspan,y0,yp0,options) also uses the integration settings defined by options, which is an argument created using the odeset function. For example, use the AbsTol and RelTol options to specify absolute and relative error tolerances, or the Mass option to provide a mass matrix.

[t,y,te,ye,ie] = ode15i(odefun,tspan,y0,yp0,options) additionally finds where functions of (t,y,y'), called event functions, are zero. In the output, te is the time of the event, ye is the solution at the time of the event, and ie is the index of the triggered event.

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 'Events' property to a function, such as myEventFcn or @myEventFcn, and creating a corresponding function: [value,isterminal,direction] = myEventFcn(t,y,yp). For more information, see ODE Event Location.

sol = ode15i(___) returns a structure that you can use with deval to evaluate the solution at any point on the interval [t0 tf]. You can use any of the input argument combinations in previous syntaxes.


collapse all

Solve Weissinger Implicit ODE

Use decic to compute consistent initial conditions for the Weissinger implicit ODE. decic holds fixed the initial value for y(t0) and computes a consistent initial value for y'(t0). The weissinger function evaluates the residual of the implicit ODE.

t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0);

Use the result returned by decic with ode15i to solve the ODE. Plot the numerical solution, y, against the analytical solution, ytrue.

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);
ytrue = sqrt(t.^2 + 0.5);

Solve Robertson Problem as Implicit Differential Algebraic Equations (DAEs)

This example reformulates a system of ODEs as a fully implicit system of differential algebraic equations (DAEs). The Robertson problem coded by hb1ode.m is a classic test problem for programs that solve stiff ODEs. The system of equations is

$$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 4y_2y_3-(3 \times 10^7)y_2^2\\ y'_3 &= (3 \times

hb1ode solves this system of ODEs to steady state with the initial conditions $y_1 = 1$, $y_2 = 0$, and $y_3 = 0$. But the equations also satisfy a linear conservation law,

$$y'_1 + y'_2 + y'_3 = 0.$$

In terms of the solution and initial conditions, the conservation law is

$$y_1 + y_2 + y_3 = 1.$$

The problem can be rewritten as a system of DAEs by using the conservation law to determine the state of $y_3$. This reformulates the problem as the implicit DAE system

$$\begin{array}{cl} 0 &= y'_1 +0.04y_1 - 10^4 y_2y_3\\ 0 &= y'_2 -0.04y_1
+ 10^4 y_2y_3+(3 \times 10^7)y_2^2\\ 0 &= y_1 + y_2 + y_3 -

The function robertsidae encodes this DAE system.

function res = robertsidae(t,y,yp)
res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);
   yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;
   y(1) + y(2) + y(3) - 1];

The full example code for this formulation of the Robertson problem is available in ihb1dae.m.

Set the error tolerances and the value of $\partial f / \partial y'$.

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...
   'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});

Use decic to compute consistent initial conditions from guesses. Fix the first two components of y0 to get the same consistent initial conditions as found by ode15s in hb1dae.m, which formulates this problem as a semi-explicit DAE system.

y0 = [1; 0; 1e-3];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);

Solve the system of DAEs using ode15i.

tspan = [0 4*logspace(-6,6)];
[t,y] = ode15i(@robertsidae,tspan,y0,yp0,options);

Plot the solution components. Since the second solution component is small relative to the others, multiply it by 1e4 before plotting.

y(:,2) = 1e4*y(:,2);
ylabel('1e4 * y(:,2)')
title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

Input Arguments

collapse all

odefun — Functions to solvefunction handle

Functions to solve, specified as a function handle that defines the functions to be integrated.

The function f = odefun(t,y,yp), for a scalar t and column vectors y and yp, must return a column vector f of data type single or double that corresponds to f(t,y,y'). odefun must accept the three inputs for t, y, and yp even if one of the inputs is not used in the function.

For example, to solve y'y=0, use this function.

function f = odefun(t,y,yp)
f = yp - y;

For a system of equations, the output of odefun is a vector. Each equation becomes an element in the solution vector. For example, to solve


use this function.

function dy = odefun(t,y,yp)
dy = zeros(2,1);
dy(1) = yp(1)-y(2);
dy(2) = yp(2)+1;

For information on how to provide additional parameters to the function odefun, see Parameterizing Functions.

Example: @myFcn

Data Types: function_handle

tspan — Interval of integrationvector

Interval of integration, specified as a vector. At minimum, tspan must be a two element vector [t0 tf] specifying the initial and final times. To obtain solutions at specific times between t0 and tf, use a longer vector of the form [t0,t1,t2,...,tf]. The elements in tspan must be all increasing or all decreasing.

The solver imposes the initial conditions, y0, at tspan(1), and then integrates from tspan(1) to tspan(end):

  • If tspan has two elements, [t0 tf], then the solver returns the solution evaluated at each internal integration step within the interval.

  • If tspan contains more than two elements [t0,t1,t2,...,tf], then the solver returns the solution evaluated at the given points. This does not affect the internal steps that the solver uses to traverse from tspan(1) to tspan(end). Thus, the solver does not necessarily step precisely to each point specified in tspan. However, the solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step.

    Specifying several intermediate points has little effect on the efficiency of computation, but for large systems it can affect memory management.

The solution obtained by the solver might be different depending on whether you specify tspan as a two-element vector or as a vector with intermediate points. If tspan contains several intermediate points, then they give an indication of the scale for the problem, which can affect the size of the initial step taken by the solver.

Example: [1 10]

Example: [1 3 5 7 9 10]

Data Types: single | double

y0 — Initial conditions for yvector

Initial conditions for y, specified as a vector. y0 must be the same length as the vector output of odefun, so that y0 contains an initial condition for each equation defined in odefun.

The initial conditions for y0 and yp0 must be consistent, meaning that f(t0,y0,y'0)=0. Use the decic function to compute consistent initial conditions close to guessed values.

Data Types: single | double

yp0 — Initial conditions for y'vector

Initial conditions for y', specified as a vector. yp0 must be the same length as the vector output of odefun, so that yp0 contains an initial condition for each variable defined in odefun.

The initial conditions for y0 and yp0 must be consistent, meaning that f(t0,y0,y'0)=0. Use the decic function to compute consistent initial conditions close to guessed values.

Data Types: single | double

options — Option structurestructure array

Option structure, specified as a structure array. Use the odeset function to create or modify the option structure.

See Summary of ODE Options for a list of which options are compatible with each ODE solver.

Example: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) specifies a relative error tolerance of 1e-5, turns on the display of solver statistics, and specifies the output function @odeplot to plot the solution as it is computed.

Data Types: struct

Output Arguments

collapse all

t — Evaluation pointscolumn vector

Evaluation points, returned as a column vector.

  • If tspan contains two elements, [t0 tf], then t contains the internal evaluation points used to perform the integration.

  • If tspan contains more than two elements, then t is the same as tspan.

y — Solutionsarray

Solutions, returned as an array. Each row in y corresponds to the solution at the value returned in the corresponding row of t.

te — Time of eventscolumn vector

Time of events, returned as a column vector. The event times in te correspond to the solutions returned in ye, and ie specifies which event occurred.

ye — Solution at time of eventsarray

Solution at time of events, returned as an array. The event times in te correspond to the solutions returned in ye, and ie specifies which event occurred.

ie — Index of vanishing event functioncolumn vector

Index of vanishing event function, returned as a column vector. The event times in te correspond to the solutions returned in ye, and ie specifies which event occurred.

sol — Structure for evaluationstructure array

Structure for evaluation, returned as a structure array. Use this structure with the deval function to evaluate the solution at any point in the interval [t0 tf]. The sol structure array always includes these fields:

Structure FieldDescription


Row vector of the steps chosen by the solver.


Solutions. Each column sol.y(:,i) contains the solution at time sol.x(i).


Solver name.

Additionally, if you specify the Events option and events are detected, then sol also includes these fields:

Structure FieldDescription


Points when events occurred. sol.xe(end) contains the exact point of a terminal event, if any.

Solutions that correspond to events in sol.xe.

Indices into the vector returned by the function specified in the Events option. The values indicate which event the solver detected.

More About

collapse all


  • Providing the Jacobian matrix to ode15i is critical for reliability and efficiency. Alternatively, if the system is large and sparse, then providing the Jacobian sparsity pattern also assists the solver. In either case, use odeset to pass in the matrices using the Jacobian or JPattern options.


ode15i is a variable-step, variable-order (VSVO) solver based on the backward differentiation formulas (BDFs) of orders 1 to 5. ode15i is designed to be used with fully implicit differential equations and index-1 differential algebraic equations (DAEs). The helper function decic computes consistent initial conditions that are suitable to be used with ode15i [1].


[1] Lawrence F. Shampine, "Solving 0 = F(t, y(t), y′(t)) in MATLAB," Journal of Numerical Mathematics, Vol.10, No.4, 2002, pp. 291-310.

See Also

| | | | |

Introduced before R2006a

Was this topic helpful?