This function solves implicit differential equations (IDE).
  The solution technique is based on the 3-stage Lobatto IIIC implicit Runge-Kutta
{Code, Test problems}

     [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)

  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)
                 |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.
           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
              is the absolute error tolerance. (Default 1.0e-6).
              is the relative error tolerance. (Default 1.0e-6).
           Note that ATOL and RTOL can also be (n x 1) vectors.
              is the initial step size. (Default 1.0e-8).
              is maximum number of steps that is performed by the function.
              (Default 1000).
              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).
              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.
  (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;
     u = 0;
  res(1) = yp(1) - y(2);
  res(2) = yp(2) + 0.1 * y(2) + 3.0 * y(1) - u;
  (2) The file idejac.m contains
  function [J, M] = idejac(y, yp, t)
  J = [0, -1; 0.1, 3.0];
  M = eye(2);
  (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);
  See also, ide_options, ide_set_option