l3cide:

This function solves implicit differential equations (IDE).
The solution technique is based on the 3-stage Lobatto IIIC implicit Runge-Kutta
coefficients.

{Code,   Usage:

[TOUT,YOUT,INFO] = l3cide('IDEFUN', 'IDEJAC', TSPAN, Y0, YP0)  or
[TOUT,YOUT,INFO] = l3cide('IDEFUN', 'IDEJAC', TSPAN, Y0, YP0, OPTIONS) or
[TOUT,YOUT,INFO] = l3cide('IDEFUN', '', TSPAN, Y0, YP0)  or
[TOUT,YOUT,INFO] = l3cide('IDEFUN', '', TSPAN, Y0, YP0, OPTIONS)

Here,

IDEFUN   is a function that returns the (n x 1) vector res = phi(y, yp, t),
where y is the (n x 1) state vector of the IDE, yp = dy/dt is the
(n x 1) state derivative, t is the time, and phi is the (n x 1)
system of IDE. This function has the form

function res = IDEFUN(y, yp, t)

IDEJAC   is a function that computes the (n x n) Jacobian J = phi/dy, and the
(nxn) matrix M = dphi/dyp.  This function has the form

function [J, M] = IDEJAC(y, yp, t)

where

|dphi_1/dy_1 dphi_1/dy_2 ... dphi_1/dy_n|
|dphi_2/dy_1 dphi_2/dy_2 ... dphi_2/dy_n|
J =  |   .                                   | and
|   .                                   |
|dphi_n/dy_1 dphi_n/dy_2 ... dphi_n/dy_n|

|dphi_1/dyp_1 dphi_1/dyp_2 ... dphi_1/dyp_n|
|dphi_2/dyp_1 dphi_2/dyp_2 ... dphi_2/dyp_n|
M =  |   .                                      |
|   .                                      |
|dphi_n/dyp_1 dphi_n/dyp_2 ... dphi_n/dyp_n|

If IDEJAC is not included in the argument list for l3cide then J and
M are computed via a finite difference approximation.

TSPAN    This vector defines the limits of integration, as well as
intermidiate time points where the solution to the IDE is computed.
Specifically,

TSPAN is a (Ntspan x 1) vector

TSPAN = [t_1; t_2; t_3; ...; t_Ntspan]

where the times t_1, t_2, ..., t_Ntspan is a monotone increasing
(or decreasing) sequence.  The limits of integration are t_1 (the
initial time, and T_Ntspan (the final time).  Solutions to the IDE
are computed a the points t_2, t_3, ..., t_Ntspan, and stored in the
matrix YOUT.  On successful termination of the function TOUT = TSPAN.

Y0       This (n x 1) vector gives the initial state for the IDE.

YP0      This (n x 1) vector gives the initial state derivative for the IDE.
Moreover, it is assumed that these initial conditions are consistent,
i.e., res = phi(Y0, YP0, TSPAN(1)) = 0.  Note that the function
terminates if ||phi(Y0, YP0, TSPAN(1))|| > sqrt(eps), where eps is
the machine precision.

OPTIONS  This is a structure that provides various options for the numerical
solution algorithm.  The default values for the members of this
structure can be ascertained using the function

OPTIONS = ide_options()

The function options can be assigned using the function

OPTIONS = ide_set_option(OPTIONS, 'OPTION', value)

The options for the function l3cide and the default values are as
follows;

OPTIONS.ATOL
is the absolute error tolerance. (Default 1.0e-6).

OPTIONS.RTOL
is the relative error tolerance. (Default 1.0e-6).

Note that ATOL and RTOL can also be (n x 1) vectors.

OPTIONS.INITIAL_STEP_SIZE
is the initial step size. (Default 1.0e-8).

OPTIONS.MAX_STEPS
is maximum number of steps that is performed by the function.
(Default 1000).

OPTIONS.DIFF_INDEX
is a (n x 1) vector that indicates the differentiation index of
the i-th state variable, y(i). (Default [], i.e., all state
variables are assumed to be index 0).

OPTIONS.T_EVENT
is a (P x 1) vector of time points of the form

OPTIONS.T_EVENT = [te_1;te_1;te_2;...;te_P]

where the times te_1, te_2, ..., te_P is a monotone increasing
(or decreasing) sequence.  In the algorithm the step sizes are
selected such that the event times, T_EVENT(k), are at the end
(or start) of the integration interval.  This is useful for modeling
discontinuous time varying inputs or, discontinuous implicit
differential equations.

TOUT     is a (Ntspan x 1) vector of output times at which the solution is
computed. (Ntspan is defined below.)

YOUT     is a (Ntspan x n) matrix of solutions to the IDE, where n is the
dimension of the IDE.  (Ntspan is defined below.)

INFO     is a structure that provides some statistics related to the solution
algorithm.  In particular,

INFO.nfun    - is the number of function evaluations.
INFO.njac    - is the number of Jacobian evaluations.
INFO.naccept - is the number of successful steps.
INFO.nreject - is the number of failed steps.

Example:

(1) The file ide.m contains

function res = ide(y, yp, t)

res = zeros(2, 1);

if ((t >= 0) && (t < 1))
u = 1;
elseif ((t >= 1) && (t < 2))
u = -1;
elseif((t >= 2) && (t < 3))
u = 1;
else
u = 0;
end;
res(1) = yp(1) - y(2);
res(2) = yp(2) + 0.1 * y(2) + 3.0 * y(1) - u;
return;

(2) The file idejac.m contains

function [J, M] = idejac(y, yp, t)

J = [0, -1; 0.1, 3.0];
M = eye(2);
return;

(3) The file runide.m contains

clear all;
tspan = linspace(0, 60.0, 600)';
l3cide_opt = ide_options();
l3cide_opt = ide_set_option(l3cide_opt, 'T_EVENT', [1; 2; 3]);
y0 = [0; 0];
yp0 = [0; 0];
yp0 = -ide(y0, yp0, 0);

[T, Y, info] = l3cide('ide', 'idejac', tspan, y0, yp0, l3cide_opt);