Differential Equations
**********************

Octave has built-in functions for solving ordinary differential equations (ODEs), and differential-algebraic equations (DAEs).

Ordinary Differential Equations
===============================

The function lsode can be used to solve ODEs of the form

     dx
     -- = f (x, t)
     dt

using Hindmarsh's ODE solver LSODE.

 -- [X, ISTATE, MSG] =: lsode (FCN, X_0, T)
 -- [X, ISTATE, MSG] =: lsode (FCN, X_0, T, T_CRIT)
     Ordinary Differential Equation (ODE) solver.

     The set of differential equations to solve is

          dx
          -- = f (x, t)
          dt

     with

          x(t_0) = x_0

     The solution is returned in the matrix X, with each row corresponding to an element of the vector T.  The first element of T should be t_0 and should correspond to the initial state of the system X_0, so that the first row of the output is X_0.

     The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of right hand sides for the set of equations.  The function must have the form

          XDOT = f (X, T)

     in which XDOT and X are vectors and T is a scalar.

     If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the Jacobian of f.  The Jacobian function must have the form

          JAC = j (X, T)

     in which JAC is the matrix of partial derivatives

                       | df_1  df_1       df_1 |
                       | ----  ----  ...  ---- |
                       | dx_1  dx_2       dx_N |
                       |                       |
                       | df_2  df_2       df_2 |
                       | ----  ----  ...  ---- |
                df_i   | dx_1  dx_2       dx_N |
          jac = ---- = |                       |
                dx_j   |  .    .     .    .    |
                       |  .    .      .   .    |
                       |  .    .       .  .    |
                       |                       |
                       | df_M  df_M       df_M |
                       | ----  ----  ...  ---- |
                       | dx_1  dx_2       dx_N |

     The second argument specifies the initial state of the system x_0.  The third argument is a vector, T, specifying the time values for which a solution is sought.

     The fourth argument is optional, and may be used to specify a set of times that the ODE solver should not integrate past.  It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative.

     After a successful computation, the value of ISTATE will be 2 (consistent with the Fortran version of LSODE).

     If the computation is not successful, ISTATE will be something other than 2 and MSG will contain additional information.

     You can use the function lsode_options to set optional parameters for lsode.

     See Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers, in Scientific Computing, R. S. Stepleman, editor, (1983) or https://computing.llnl.gov/projects/odepack for more information about the inner workings of lsode.

     Example: Solve the Van der Pol equation

          fvdp = @(Y,T) [Y(2); (1 - Y(1)^2) * Y(2) - Y(1)];
          T = linspace (0, 20, 100);
          Y = lsode (fvdp, [2; 0], T);

     daspk, dassl, dasrt

 -- lsode_options: ()
 -- val =: lsode_options (OPT)
 -- lsode_options: (OPT, VAL)
     Query or set options for the function lsode.

     When called with no arguments, the names of all available options and their current values are displayed.

     Given one argument, return the value of the option OPT.

     When called with two arguments, lsode_options sets the option OPT to value VAL.

     Options include

     "absolute tolerance"
          Absolute tolerance.  May be either vector or scalar.  If a vector, it must match the dimension of the state vector.

     "relative tolerance"
          Relative tolerance parameter.  Unlike the absolute tolerance, this parameter may only be a scalar.

          The local error test applied at each integration step is

                 abs (local error in x(i)) <= ...
                     rtol * abs (y(i)) + atol(i)

     "integration method"
          A string specifying the method of integration to use to solve the ODE system.  Valid values are

          "adams"
          "non-stiff"
               No Jacobian used (even if it is available).

          "bdf"
          "stiff"
               Use stiff backward differentiation formula (BDF) method.  If a function to compute the Jacobian is not supplied, lsode will compute a finite difference approximation of the Jacobian matrix.

     "initial step size"
          The step size to be attempted on the first step (default is determined automatically).

     "maximum order"
          Restrict the maximum order of the solution method.  If using the Adams method, this option must be between 1 and 12.  Otherwise, it must be between 1 and 5, inclusive.

     "maximum step size"
          Setting the maximum stepsize will avoid passing over very large regions  (default is not specified).

     "minimum step size"
          The minimum absolute step size allowed (default is 0).

     "step limit"
          Maximum number of steps allowed (default is 100000).

     "jacobian type"
          A string specifying the type of Jacobian used with the stiff backward differentiation formula (BDF) integration method.  Valid values are

          "full"
               The default.  All partial derivatives are approximated or used from the user-supplied Jacobian function.

          "banded"
               Only the diagonal and the number of lower and upper subdiagonals specified by the options "lower jacobian subdiagonals" and "upper jacobian subdiagonals", respectively, are approximated or used from the user-supplied Jacobian function.  A user-supplied Jacobian function may set all other partial derivatives to arbitrary values.

          "diagonal"
               If a Jacobian function is supplied by the user, this setting has no effect. A Jacobian approximated by lsode is restricted to the diagonal, where each partial derivative is computed by applying a finite change to all elements of the state together; if the real Jacobian is indeed always diagonal, this has the same effect as applying the finite change only to the respective element of the state, but is more efficient.

     "lower jacobian subdiagonals"
          Number of lower subdiagonals used if option "jacobian type" is set to "banded".  The default is zero.

     "upper jacobian subdiagonals"
          Number of upper subdiagonals used if option "jacobian type" is set to "banded".  The default is zero.

Here is an example of solving a set of three differential equations using lsode.  Given the function

     ## oregonator differential equation
     function xdot = f (x, t)

       xdot = zeros (3,1);

       xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) ...
                 - 8.375e-06*x(1)^2);
       xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
       xdot(3) = 0.161*(x(1) - x(3));

     endfunction

and the initial condition x0 = [ 4; 1.1; 4 ], the set of equations can be integrated using the command

     t = linspace (0, 500, 1000);

     y = lsode ("f", x0, t);

If you try this, you will see that the value of the result changes dramatically between T = 0 and 5, and again around T = 305. A more efficient set of output points might be

     t = [0, logspace(-1, log10(303), 150), ...
             logspace(log10(304), log10(500), 150)];

An m-file for the differential equation used above is included with the Octave distribution in the examples directory under the name oregonator.m.

Differential-Algebraic Equations
================================

The function daspk can be used to solve DAEs of the form

     0 = f (x-dot, x, t),    x(t=0) = x_0, x-dot(t=0) = x-dot_0

where x-dot is the derivative of x.  The equation is solved using Petzold's DAE solver DASPK.

 -- [X, XDOT, ISTATE, MSG] =: daspk (FCN, X_0, XDOT_0, T, T_CRIT)
     Solve a set of differential-algebraic equations.

     daspk solves the set of equations

          0 = f (x, xdot, t)

     with

          x(t_0) = x_0, xdot(t_0) = xdot_0

     The solution is returned in the matrices X and XDOT, with each row in the result matrices corresponding to one of the elements in the vector T.  The first element of T should be t_0 and correspond to the initial state of the system X_0 and its derivative XDOT_0, so that the first row of the output X is X_0 and the first row of the output XDOT is XDOT_0.

     The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of residuals for the set of equations.  It must have the form

          RES = f (X, XDOT, T)

     in which X, XDOT, and RES are vectors, and T is a scalar.

     If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the modified Jacobian

                df       df
          jac = -- + c ------
                dx     d xdot

     The modified Jacobian function must have the form

          JAC = j (X, XDOT, T, C)

     The second and third arguments to daspk specify the initial condition of the states and their derivatives, and the fourth argument specifies a vector of output times at which the solution is desired, including the time corresponding to the initial condition.

     The set of initial states and derivatives are not strictly required to be consistent.  If they are not consistent, you must use the daspk_options function to provide additional information so that daspk can compute a consistent starting point.

     The fifth argument is optional, and may be used to specify a set of times that the DAE solver should not integrate past.  It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative.

     After a successful computation, the value of ISTATE will be greater than zero (consistent with the Fortran version of DASPK).

     If the computation is not successful, the value of ISTATE will be less than zero and MSG will contain additional information.

     You can use the function daspk_options to set optional parameters for daspk. dassl

 -- daspk_options: ()
 -- val =: daspk_options (OPT)
 -- daspk_options: (OPT, VAL)
     Query or set options for the function daspk.

     When called with no arguments, the names of all available options and their current values are displayed.

     Given one argument, return the value of the option OPT.

     When called with two arguments, daspk_options sets the option OPT to value VAL.

     Options include

     "absolute tolerance"
          Absolute tolerance.  May be either vector or scalar.  If a vector, it must match the dimension of the state vector, and the relative tolerance must also be a vector of the same length.

     "relative tolerance"
          Relative tolerance.  May be either vector or scalar.  If a vector, it must match the dimension of the state vector, and the absolute tolerance must also be a vector of the same length.

          The local error test applied at each integration step is

                 abs (local error in x(i))
                      <= rtol(i) * abs (Y(i)) + atol(i)

     "compute consistent initial condition"
          Denoting the differential variables in the state vector by Y_d and the algebraic variables by Y_a, ddaspk can solve one of two initialization problems:

           1. Given Y_d, calculate Y_a and Y'_d

           2. Given Y', calculate Y.

          In either case, initial values for the given components are input, and initial guesses for the unknown components must also be provided as input.  Set this option to 1 to solve the first problem, or 2 to solve the second (the default is 0, so you must provide a set of initial conditions that are consistent).

          If this option is set to a nonzero value, you must also set the "algebraic variables" option to declare which variables in the problem are algebraic.

     "use initial condition heuristics"
          Set to a nonzero value to use the initial condition heuristics options described below.

     "initial condition heuristics"
          A vector of the following parameters that can be used to control the initial condition calculation.

          MXNIT
               Maximum number of Newton iterations (default is 5).

          MXNJ
               Maximum number of Jacobian evaluations (default is 6).

          MXNH
               Maximum number of values of the artificial stepsize parameter to be tried if the "compute consistent initial condition" option has been set to 1 (default is 5).

               Note that the maximum total number of Newton iterations allowed is MXNIT*MXNJ*MXNH if the "compute consistent initial condition" option has been set to 1 and MXNIT*MXNJ if it is set to 2.

          LSOFF
               Set to a nonzero value to disable the linesearch algorithm (default is 0).

          STPTOL
               Minimum scaled step in linesearch algorithm (default is eps^(2/3)).

          EPINIT
               Swing factor in the Newton iteration convergence test.  The test is applied to the residual vector, premultiplied by the approximate Jacobian.  For convergence, the weighted RMS norm of this vector (scaled by the error weights) must be less than EPINIT*EPCON, where EPCON = 0.33 is the analogous test constant used in the time steps.  The default is EPINIT = 0.01.

     "print initial condition info"
          Set this option to a nonzero value to display detailed information about the initial condition calculation (default is 0).

     "exclude algebraic variables from error test"
          Set to a nonzero value to exclude algebraic variables from the error test.  You must also set the "algebraic variables" option to declare which variables in the problem are algebraic (default is 0).

     "algebraic variables"
          A vector of the same length as the state vector.  A nonzero element indicates that the corresponding element of the state vector is an algebraic variable (i.e., its derivative does not appear explicitly in the equation set).

          This option is required by the "compute consistent initial condition" and "exclude algebraic variables from error test" options.

     "enforce inequality constraints"
          Set to one of the following values to enforce the inequality constraints specified by the "inequality constraint types" option (default is 0).

           1. To have constraint checking only in the initial condition calculation.

           2. To enforce constraint checking during the integration.

           3. To enforce both options 1 and 2.

     "inequality constraint types"
          A vector of the same length as the state specifying the type of inequality constraint.  Each element of the vector corresponds to an element of the state and should be assigned one of the following codes

          -2
               Less than zero.

          -1
               Less than or equal to zero.

          0
               Not constrained.

          1
               Greater than or equal to zero.

          2
               Greater than zero.

          This option only has an effect if the "enforce inequality constraints" option is nonzero.

     "initial step size"
          Differential-algebraic problems may occasionally suffer from severe scaling difficulties on the first step.  If you know a great deal about the scaling of your problem, you can help to alleviate this problem by specifying an initial stepsize (default is computed automatically).

     "maximum order"
          Restrict the maximum order of the solution method.  This option must be between 1 and 5, inclusive (default is 5).

     "maximum step size"
          Setting the maximum stepsize will avoid passing over very large regions (default is not specified).

Octave also includes DASSL, an earlier version of DASPK, and DASRT, which can be used to solve DAEs with constraints (stopping conditions).

 -- [X, XDOT, ISTATE, MSG] =: dassl (FCN, X_0, XDOT_0, T, T_CRIT)
     Solve a set of differential-algebraic equations.

     dassl solves the set of equations

          0 = f (x, xdot, t)

     with

          x(t_0) = x_0, xdot(t_0) = xdot_0

     The solution is returned in the matrices X and XDOT, with each row in the result matrices corresponding to one of the elements in the vector T.  The first element of T should be t_0 and correspond to the initial state of the system X_0 and its derivative XDOT_0, so that the first row of the output X is X_0 and the first row of the output XDOT is XDOT_0.

     The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of residuals for the set of equations.  It must have the form

          RES = f (X, XDOT, T)

     in which X, XDOT, and RES are vectors, and T is a scalar.

     If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the modified Jacobian

                df       df
          jac = -- + c ------
                dx     d xdot

     The modified Jacobian function must have the form

          JAC = j (X, XDOT, T, C)

     The second and third arguments to dassl specify the initial condition of the states and their derivatives, and the fourth argument specifies a vector of output times at which the solution is desired, including the time corresponding to the initial condition.

     The set of initial states and derivatives are not strictly required to be consistent.  In practice, however, DASSL is not very good at determining a consistent set for you, so it is best if you ensure that the initial values result in the function evaluating to zero.

     The fifth argument is optional, and may be used to specify a set of times that the DAE solver should not integrate past.  It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative.

     After a successful computation, the value of ISTATE will be greater than zero (consistent with the Fortran version of DASSL).

     If the computation is not successful, the value of ISTATE will be less than zero and MSG will contain additional information.

     You can use the function dassl_options to set optional parameters for dassl. daspk, dasrt, lsode

 -- dassl_options: ()
 -- val =: dassl_options (OPT)
 -- dassl_options: (OPT, VAL)
     Query or set options for the function dassl.

     When called with no arguments, the names of all available options and their current values are displayed.

     Given one argument, return the value of the option OPT.

     When called with two arguments, dassl_options sets the option OPT to value VAL.

     Options include

     "absolute tolerance"
          Absolute tolerance.  May be either vector or scalar.  If a vector, it must match the dimension of the state vector, and the relative tolerance must also be a vector of the same length.

     "relative tolerance"
          Relative tolerance.  May be either vector or scalar.  If a vector, it must match the dimension of the state vector, and the absolute tolerance must also be a vector of the same length.

          The local error test applied at each integration step is

                 abs (local error in x(i))
                      <= rtol(i) * abs (Y(i)) + atol(i)

     "compute consistent initial condition"
          If nonzero, dassl will attempt to compute a consistent set of initial conditions.  This is generally not reliable, so it is best to provide a consistent set and leave this option set to zero.

     "enforce nonnegativity constraints"
          If you know that the solutions to your equations will always be non-negative, it may help to set this parameter to a nonzero value.  However, it is probably best to try leaving this option set to zero first, and only setting it to a nonzero value if that doesn't work very well.

     "initial step size"
          Differential-algebraic problems may occasionally suffer from severe scaling difficulties on the first step.  If you know a great deal about the scaling of your problem, you can help to alleviate this problem by specifying an initial stepsize.

     "maximum order"
          Restrict the maximum order of the solution method.  This option must be between 1 and 5, inclusive.

     "maximum step size"
          Setting the maximum stepsize will avoid passing over very large regions  (default is not specified).

     "step limit"
          Maximum number of integration steps to attempt on a single call to the underlying Fortran code.

 -- [X, XDOT, T_OUT, ISTAT, MSG] =: dasrt (FCN, G, X_0, XDOT_0, T)
 -- ... =: dasrt (FCN, G, X_0, XDOT_0, T, T_CRIT)
 -- ... =: dasrt (FCN, X_0, XDOT_0, T)
 -- ... =: dasrt (FCN, X_0, XDOT_0, T, T_CRIT)
     Solve a set of differential-algebraic equations.

     dasrt solves the set of equations

          0 = f (x, xdot, t)

     with

          x(t_0) = x_0, xdot(t_0) = xdot_0

     with functional stopping criteria (root solving).

     The solution is returned in the matrices X and XDOT, with each row in the result matrices corresponding to one of the elements in the vector T_OUT.  The first element of T should be t_0 and correspond to the initial state of the system X_0 and its derivative XDOT_0, so that the first row of the output X is X_0 and the first row of the output XDOT is XDOT_0.

     The vector T provides an upper limit on the length of the integration.  If the stopping condition is met, the vector T_OUT will be shorter than T, and the final element of T_OUT will be the point at which the stopping condition was met, and may not correspond to any element of the vector T.

     The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of residuals for the set of equations.  It must have the form

          RES = f (X, XDOT, T)

     in which X, XDOT, and RES are vectors, and T is a scalar.

     If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the modified Jacobian

                df       df
          jac = -- + c ------
                dx     d xdot

     The modified Jacobian function must have the form

          JAC = j (X, XDOT, T, C)

     The optional second argument names a function that defines the constraint functions whose roots are desired during the integration. This function must have the form

          G_OUT = g (X, T)

     and return a vector of the constraint function values. If the value of any of the constraint functions changes sign, DASRT will attempt to stop the integration at the point of the sign change.

     If the name of the constraint function is omitted, dasrt solves the same problem as daspk or dassl.

     Note that because of numerical errors in the constraint functions due to round-off and integration error, DASRT may return false roots, or return the same root at two or more nearly equal values of T.  If such false roots are suspected, the user should consider smaller error tolerances or higher precision in the evaluation of the constraint functions.

     If a root of some constraint function defines the end of the problem, the input to DASRT should nevertheless allow integration to a point slightly past that root, so that DASRT can locate the root by interpolation.

     The third and fourth arguments to dasrt specify the initial condition of the states and their derivatives, and the fourth argument specifies a vector of output times at which the solution is desired, including the time corresponding to the initial condition.

     The set of initial states and derivatives are not strictly required to be consistent.  In practice, however, DASSL is not very good at determining a consistent set for you, so it is best if you ensure that the initial values result in the function evaluating to zero.

     The sixth argument is optional, and may be used to specify a set of times that the DAE solver should not integrate past.  It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative.

     After a successful computation, the value of ISTATE will be greater than zero (consistent with the Fortran version of DASSL).

     If the computation is not successful, the value of ISTATE will be less than zero and MSG will contain additional information.

     You can use the function dasrt_options to set optional parameters for dasrt. dasrt_options, daspk, dasrt, lsode

 -- dasrt_options: ()
 -- val =: dasrt_options (OPT)
 -- dasrt_options: (OPT, VAL)
     Query or set options for the function dasrt.

     When called with no arguments, the names of all available options and their current values are displayed.

     Given one argument, return the value of the option OPT.

     When called with two arguments, dasrt_options sets the option OPT to value VAL.

     Options include

     "absolute tolerance"
          Absolute tolerance.  May be either vector or scalar.  If a vector, it must match the dimension of the state vector, and the relative tolerance must also be a vector of the same length.

     "relative tolerance"
          Relative tolerance.  May be either vector or scalar.  If a vector, it must match the dimension of the state vector, and the absolute tolerance must also be a vector of the same length.

          The local error test applied at each integration step is

                 abs (local error in x(i)) <= ...
                     rtol(i) * abs (Y(i)) + atol(i)

     "initial step size"
          Differential-algebraic problems may occasionally suffer from severe scaling difficulties on the first step.  If you know a great deal about the scaling of your problem, you can help to alleviate this problem by specifying an initial stepsize.

     "maximum order"
          Restrict the maximum order of the solution method.  This option must be between 1 and 5, inclusive.

     "maximum step size"
          Setting the maximum stepsize will avoid passing over very large regions.

     "step limit"
          Maximum number of integration steps to attempt on a single call to the underlying Fortran code.

See K. E. Brenan, et al., Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, North-Holland (1989), DOI: https://doi.org/10.1137/1.9781611971224, for more information about the implementation of DASSL.

Matlab-compatible solvers
=========================

Octave also provides a set of solvers for initial value problems for ordinary differential equations (ODEs) that have a MATLAB-compatible interface. The options for this class of methods are set using the functions.

  * odeset

  * odeget

Currently implemented solvers are:

  * Runge-Kutta methods

       * ode45 integrates a system of non-stiff ODEs or index-1 differential-algebraic equations (DAEs) using the high-order, variable-step Dormand-Prince method.  It requires six function evaluations per integration step, but may take larger steps on smooth problems than ode23: potentially offering improved efficiency at smaller tolerances.

       * ode23 integrates a system of non-stiff ODEs or (or index-1 DAEs).  It uses the third-order Bogacki-Shampine method and adapts the local step size in order to satisfy a user-specified tolerance.  The solver requires three function evaluations per integration step.

       * ode23s integrates a system of stiff ODEs (or index-1 DAEs) using a modified second-order Rosenbrock method.

  * Linear multistep methods

       * ode15s integrates a system of stiff ODEs (or index-1 DAEs) using a variable step, variable order method based on Backward Difference Formulas (BDF).

       * ode15i integrates a system of fully-implicit ODEs (or index-1 DAEs) using the same variable step, variable order method as ode15s.  decic can be used to compute consistent initial conditions for ode15i.

Detailed information on the solvers are given in L. F. Shampine and M. W. Reichelt, The MATLAB ODE Suite, SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22, DOI: https://doi.org/10.1137/S1064827594276424.

 -- [T, Y] =: ode45 (FCN, TRANGE, INIT)
 -- [T, Y] =: ode45 (FCN, TRANGE, INIT, ODE_OPT)
 -- [T, Y, TE, YE, IE] =: ode45 (...)
 -- SOLUTION =: ode45 (...)
 -- ode45: (...)

     Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Dormand-Prince method of order 4.

     FCN is a function handle, inline function, or string containing the name of the function that defines the ODE: y' = f(t,y).  The function must accept two inputs where the first is time T and the second is a column vector of unknowns Y.

     TRANGE specifies the time interval over which the ODE will be evaluated.  Typically, it is a two-element vector specifying the initial and final times ([tinit, tfinal]).  If there are more than two elements then the solution will also be evaluated at these intermediate time instances.

     By default, ode45 uses an adaptive timestep with the integrate_adaptive algorithm.  The tolerance for the timestep computation may be changed by using the options "RelTol" and "AbsTol".

     INIT contains the initial value for the unknowns.  If it is a row vector then the solution Y will be a matrix in which each column is the solution for the corresponding initial value in INIT.

     The optional fourth argument ODE_OPT specifies non-default options to the ODE solver.  It is a structure generated by odeset.

     The function typically returns two outputs.  Variable T is a column vector and contains the times where the solution was found.  The output Y is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in T.

     The output can also be returned as a structure SOLUTION which has a field X containing a row vector of times where the solution was evaluated and a field Y containing the solution matrix such that each column corresponds to a time in X.  Use fieldnames (SOLUTION) to see the other fields and additional information returned.

     If no output arguments are requested, and no "OutputFcn" is specified in ODE_OPT, then the "OutputFcn" is set to odeplot and the results of the solver are plotted immediately.

     If using the "Events" option then three additional outputs may be returned.  TE holds the time when an Event function returned a zero. YE holds the value of the solution at time TE.  IE contains an index indicating which Event function was triggered in the case of multiple Event functions.

     Example: Solve the Van der Pol equation

          fvdp = @(T,Y) [Y(2); (1 - Y(1)^2) * Y(2) - Y(1)];
          [T,Y] = ode45 (fvdp, [0, 20], [2, 0]);

     odeset, odeget, ode23, ode15s

 -- [T, Y] =: ode23 (FCN, TRANGE, INIT)
 -- [T, Y] =: ode23 (FCN, TRANGE, INIT, ODE_OPT)
 -- [T, Y, TE, YE, IE] =: ode23 (...)
 -- SOLUTION =: ode23 (...)
 -- ode23: (...)

     Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Bogacki-Shampine method of order 3.

     FCN is a function handle, inline function, or string containing the name of the function that defines the ODE: y' = f(t,y).  The function must accept two inputs where the first is time T and the second is a column vector of unknowns Y.

     TRANGE specifies the time interval over which the ODE will be evaluated.  Typically, it is a two-element vector specifying the initial and final times ([tinit, tfinal]).  If there are more than two elements then the solution will also be evaluated at these intermediate time instances.

     By default, ode23 uses an adaptive timestep with the integrate_adaptive algorithm.  The tolerance for the timestep computation may be changed by using the options "RelTol" and "AbsTol".

     INIT contains the initial value for the unknowns.  If it is a row vector then the solution Y will be a matrix in which each column is the solution for the corresponding initial value in INIT.

     The optional fourth argument ODE_OPT specifies non-default options to the ODE solver.  It is a structure generated by odeset.

     The function typically returns two outputs.  Variable T is a column vector and contains the times where the solution was found.  The output Y is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in T.

     The output can also be returned as a structure SOLUTION which has a field X containing a row vector of times where the solution was evaluated and a field Y containing the solution matrix such that each column corresponds to a time in X.  Use fieldnames (SOLUTION) to see the other fields and additional information returned.

     If no output arguments are requested, and no "OutputFcn" is specified in ODE_OPT, then the "OutputFcn" is set to odeplot and the results of the solver are plotted immediately.

     If using the "Events" option then three additional outputs may be returned.  TE holds the time when an Event function returned a zero. YE holds the value of the solution at time TE.  IE contains an index indicating which Event function was triggered in the case of multiple Event functions.

     Example: Solve the Van der Pol equation

          fvdp = @(T,Y) [Y(2); (1 - Y(1)^2) * Y(2) - Y(1)];
          [T,Y] = ode23 (fvdp, [0, 20], [2, 0]);

     Reference: For the definition of this method see https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods. odeset, odeget, ode45, ode15s

 -- [T, Y] =: ode23s (FCN, TRANGE, INIT)
 -- [T, Y] =: ode23s (FCN, TRANGE, INIT, ODE_OPT)
 -- [T, Y] =: ode23s (..., PAR1, PAR2, ...)
 -- [T, Y, TE, YE, IE] =: ode23s (...)
 -- SOLUTION =: ode23s (...)

     Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a Rosenbrock method of order (2,3).

     FCN is a function handle, inline function, or string containing the name of the function that defines the ODE: M y' = f(t,y).  The function must accept two inputs where the first is time T and the second is a column vector of unknowns Y.  M is a constant mass matrix, non-singular and possibly sparse.  Set the field "Mass" in ODEOPTS using ODESET to specify a mass matrix.

     TRANGE specifies the time interval over which the ODE will be evaluated.  Typically, it is a two-element vector specifying the initial and final times ([tinit, tfinal]).  If there are more than two elements then the solution will also be evaluated at these intermediate time instances using an interpolation procedure of the same order as the one of the solver.

     By default, ode23s uses an adaptive timestep with the integrate_adaptive algorithm.  The tolerance for the timestep computation may be changed by using the options "RelTol" and "AbsTol".

     INIT contains the initial value for the unknowns.  If it is a row vector then the solution Y will be a matrix in which each column is the solution for the corresponding initial value in INIT.

     The optional fourth argument ODE_OPT specifies non-default options to the ODE solver.  It is a structure generated by odeset. ode23s will ignore the following options: "BDF", "InitialSlope", "MassSingular", "MStateDependence", "MvPattern", "MaxOrder", "Non-negative".

     The function typically returns two outputs.  Variable T is a column vector and contains the times where the solution was found.  The output Y is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in T.  If TRANGE specifies intermediate time steps, only those will be returned.

     The output can also be returned as a structure SOLUTION which has a field X containing a row vector of times where the solution was evaluated and a field Y containing the solution matrix such that each column corresponds to a time in X.  Use fieldnames (SOLUTION) to see the other fields and additional information returned.

     If using the "Events" option then three additional outputs may be returned.  TE holds the time when an Event function returned a zero. YE holds the value of the solution at time TE.  IE contains an index indicating which Event function was triggered in the case of multiple Event functions.

     Example: Solve the stiff Van der Pol equation

          f = @(T,Y) [Y(2); 1000*(1 - Y(1)^2) * Y(2) - Y(1)];
          opt = odeset ('Mass', [1 0; 0 1], 'MaxStep', 1e-1);
          [vt, vy] = ode23s (f, [0 2000], [2 0], opt);

     odeset, daspk, dassl

 -- [T, Y] =: ode15s (FCN, TRANGE, Y0)
 -- [T, Y] =: ode15s (FCN, TRANGE, Y0, ODE_OPT)
 -- [T, Y, TE, YE, IE] =: ode15s (...)
 -- SOLUTION =: ode15s (...)
 -- ode15s: (...)
     Solve a set of stiff Ordinary Differential Equations (ODEs) or stiff semi-explicit index 1 Differential Algebraic Equations (DAEs).

     ode15s uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.

     FCN is a function handle, inline function, or string containing the name of the function that defines the ODE: y' = f(t,y).  The function must accept two inputs where the first is time T and the second is a column vector of unknowns Y.

     TRANGE specifies the time interval over which the ODE will be evaluated.  Typically, it is a two-element vector specifying the initial and final times ([tinit, tfinal]).  If there are more than two elements then the solution will also be evaluated at these intermediate time instances.

     INIT contains the initial value for the unknowns.  If it is a row vector then the solution Y will be a matrix in which each column is the solution for the corresponding initial value in INIT.

     The optional fourth argument ODE_OPT specifies non-default options to the ODE solver.  It is a structure generated by odeset.

     The function typically returns two outputs.  Variable T is a column vector and contains the times where the solution was found.  The output Y is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in T.

     The output can also be returned as a structure SOLUTION which has a field X containing a row vector of times where the solution was evaluated and a field Y containing the solution matrix such that each column corresponds to a time in X.  Use fieldnames (SOLUTION) to see the other fields and additional information returned.

     If no output arguments are requested, and no "OutputFcn" is specified in ODE_OPT, then the "OutputFcn" is set to odeplot and the results of the solver are plotted immediately.

     If using the "Events" option then three additional outputs may be returned.  TE holds the time when an Event function returned a zero. YE holds the value of the solution at time TE.  IE contains an index indicating which Event function was triggered in the case of multiple Event functions.

     Example: Solve Robertson's equations:

          function r = robertson_dae (T, Y)
            r = [ -0.04*Y(1) + 1e4*Y(2)*Y(3)
                   0.04*Y(1) - 1e4*Y(2)*Y(3) - 3e7*Y(2)^2
          Y(1) + Y(2) + Y(3) - 1 ];
          endfunction
          opt = odeset ("Mass", [1 0 0; 0 1 0; 0 0 0], "MStateDependence", "none");
          [T,Y] = ode15s (@robertson_dae, [0, 1e3], [1; 0; 0], opt);

     decic, odeset, odeget, ode23, ode45

 -- [T, Y] =: ode15i (FCN, TRANGE, Y0, YP0)
 -- [T, Y] =: ode15i (FCN, TRANGE, Y0, YP0, ODE_OPT)
 -- [T, Y, TE, YE, IE] =: ode15i (...)
 -- SOLUTION =: ode15i (...)
 -- ode15i: (...)
     Solve a set of fully-implicit Ordinary Differential Equations (ODEs) or index 1 Differential Algebraic Equations (DAEs).

     ode15i uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.

     FCN is a function handle, inline function, or string containing the name of the function that defines the ODE: 0 = f(t,y,yp).  The function must accept three inputs where the first is time T, the second is the function value Y (a column vector), and the third is the derivative value YP (a column vector).

     TRANGE specifies the time interval over which the ODE will be evaluated.  Typically, it is a two-element vector specifying the initial and final times ([tinit, tfinal]).  If there are more than two elements then the solution will also be evaluated at these intermediate time instances.

     Y0 and YP0 contain the initial values for the unknowns Y and YP.  If they are row vectors then the solution Y will be a matrix in which each column is the solution for the corresponding initial value in Y0 and YP0.

     Y0 and YP0 must be consistent initial conditions, meaning that f(t,y0,yp0) = 0 is satisfied.  The function decic may be used to compute consistent initial conditions given initial guesses.

     The optional fifth argument ODE_OPT specifies non-default options to the ODE solver.  It is a structure generated by odeset.

     The function typically returns two outputs.  Variable T is a column vector and contains the times where the solution was found.  The output Y is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in T.

     The output can also be returned as a structure SOLUTION which has a field X containing a row vector of times where the solution was evaluated and a field Y containing the solution matrix such that each column corresponds to a time in X.  Use fieldnames (SOLUTION) to see the other fields and additional information returned.

     If no output arguments are requested, and no "OutputFcn" is specified in ODE_OPT, then the "OutputFcn" is set to odeplot and the results of the solver are plotted immediately.

     If using the "Events" option then three additional outputs may be returned.  TE holds the time when an Event function returned a zero. YE holds the value of the solution at time TE.  IE contains an index indicating which Event function was triggered in the case of multiple Event functions.

     Example: Solve Robertson's equations:

          function r = robertson_dae (T, Y, YP)
            r = [ -(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 ];
          endfunction
          [T,Y] = ode15i (@robertson_dae, [0, 1e3], [1; 0; 0], [-1e-4; 1e-4; 0]);

     decic, odeset, odeget

 -- [Y0_NEW, YP0_NEW] =: decic (FCN, T0, Y0, FIXED_Y0, YP0, FIXED_YP0)
 -- [Y0_NEW, YP0_NEW] =: decic (FCN, T0, Y0, FIXED_Y0, YP0, FIXED_YP0, OPTIONS)
 -- [Y0_NEW, YP0_NEW, RESNORM] =: decic (...)

     Compute consistent implicit ODE initial conditions Y0_NEW and YP0_NEW given initial guesses Y0 and YP0.

     A maximum of length (Y0) components between FIXED_Y0 and FIXED_YP0 may be chosen as fixed values.

     FCN is a function handle.  The function must accept three inputs where the first is time T, the second is a column vector of unknowns Y, and the third is a column vector of unknowns YP.

     T0 is the initial time such that FCN(T0, Y0_NEW, YP0_NEW) = 0, specified as a scalar.

     Y0 is a vector used as the initial guess for Y.

     FIXED_Y0 is a vector which specifies the components of Y0 to hold fixed.  Choose a maximum of length (Y0) components between FIXED_Y0 and FIXED_YP0 as fixed values. Set FIXED_Y0(i) component to 1 if you want to fix the value of Y0(i). Set FIXED_Y0(i) component to 0 if you want to allow the value of Y0(i) to change.

     YP0 is a vector used as the initial guess for YP.

     FIXED_YP0 is a vector which specifies the components of YP0 to hold fixed.  Choose a maximum of length (YP0) components between FIXED_Y0 and FIXED_YP0 as fixed values. Set FIXED_YP0(i) component to 1 if you want to fix the value of YP0(i). Set FIXED_YP0(i) component to 0 if you want to allow the value of YP0(i) to change.

     The optional seventh argument OPTIONS is a structure array.  Use odeset to generate this structure.  The relevant options are RelTol and AbsTol which specify the error thresholds used to compute the initial conditions.

     The function typically returns two outputs.  Variable Y0_NEW is a column vector and contains the consistent initial value of Y.  The output YP0_NEW is a column vector and contains the consistent initial value of YP.

     The optional third output RESNORM is the norm of the vector of residuals.  If RESNORM is small, decic has successfully computed the initial conditions.  If the value of RESNORM is large, use RelTol and AbsTol to adjust it.

     Example: Compute initial conditions for Robertson's equations:

          function r = robertson_dae (T, Y, YP)
            r = [ -(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 ];
          endfunction

          [Y0_NEW,YP0_NEW] = decic (@robertson_dae, 0, [1; 0; 0], [1; 1; 0],
          [-1e-4; 1; 0], [0; 0; 0]);

     ode15i, odeset

 -- ODESTRUCT =: odeset ()
 -- ODESTRUCT =: odeset ("FIELD1", VALUE1, "FIELD2", VALUE2, ...)
 -- ODESTRUCT =: odeset (OLDSTRUCT, "FIELD1", VALUE1, "FIELD2", VALUE2, ...)
 -- ODESTRUCT =: odeset (OLDSTRUCT, NEWSTRUCT)
 -- odeset: ()

     Create or modify an ODE options structure.

     When called with no input argument and one output argument, return a new ODE options structure that contains all possible fields initialized to their default values.  If no output argument is requested, display a list of the common ODE solver options along with their default value.

     If called with name-value input argument pairs "FIELD1", "VALUE1", "FIELD2", "VALUE2", ... return a new ODE options structure with all the most common option fields initialized, and set the values of the fields "FIELD1", "FIELD2", ... to the values VALUE1, VALUE2, ....

     If called with an input structure OLDSTRUCT then overwrite the values of the options "FIELD1", "FIELD2", ... with new values VALUE1, VALUE2, ... and return the modified structure.

     When called with two input ODE options structures OLDSTRUCT and NEWSTRUCT overwrite all values from the structure OLDSTRUCT with new values from the structure NEWSTRUCT. Empty values in NEWSTRUCT will not overwrite values in OLDSTRUCT.

     The most commonly used ODE options, which are always assigned a value by odeset, are the following:

     AbsTol: positive scalar | vector, def. 1e-6
          Absolute error tolerance.

     BDF: {"off"} | "on"
          Use BDF formulas in implicit multistep methods. Note: This option is not yet implemented.

     Events: function_handle
          Event function.  An event function must have the form [value, isterminal, direction] = my_events_f (t, y)

     InitialSlope: vector
          Consistent initial slope vector for DAE solvers.

     InitialStep: positive scalar
          Initial time step size.

     Jacobian: matrix | function_handle
          Jacobian matrix, specified as a constant matrix or a function of time and state.

     JConstant: {"off"} | "on"
          Specify whether the Jacobian is a constant matrix or depends on the state.

     JPattern: sparse matrix
          If the Jacobian matrix is sparse and non-constant but maintains a constant sparsity pattern, specify the sparsity pattern.

     Mass: matrix | function_handle
          Mass matrix, specified as a constant matrix or a function of time and state.

     MassSingular: {"maybe"} | "yes" | "on"
          Specify whether the mass matrix is singular.

     MaxOrder: {5} | 4 | 3 | 2 | 1
          Maximum order of formula.

     MaxStep: positive scalar
          Maximum time step value.

     MStateDependence: {"weak"} | "none" | "strong"
          Specify whether the mass matrix depends on the state or only on time.

     MvPattern: sparse matrix
          If the mass matrix is sparse and non-constant but maintains a constant sparsity pattern, specify the sparsity pattern. Note: This option is not yet implemented.

     NonNegative: scalar | vector
          Specify elements of the state vector that are expected to remain non-negative during the simulation.

     NormControl: {"off"} | "on"
          Control error relative to the 2-norm of the solution, rather than its absolute value.

     OutputFcn: function_handle
          Function to monitor the state during the simulation.  For the form of the function to use see odeplot.

     OutputSel: scalar | vector
          Indices of elements of the state vector to be passed to the output monitoring function.

     Refine: positive scalar
          Specify whether output should be returned only at the end of each time step or also at intermediate time instances.  The value should be a scalar indicating the number of equally spaced time points to use within each timestep at which to return output.

     RelTol: positive scalar
          Relative error tolerance.

     Stats: {"off"} | "on"
          Print solver statistics after simulation.

     Vectorized: {"off"} | "on"
          Specify whether odefcn can be passed multiple values of the state at once.

     Field names that are not in the above list are also accepted and added to the result structure.

     odeget

 -- VAL =: odeget (ODE_OPT, FIELD)
 -- VAL =: odeget (ODE_OPT, FIELD, DEFAULT)

     Query the value of the property FIELD in the ODE options structure ODE_OPT.

     If called with two input arguments and the first input argument ODE_OPT is an ODE option structure and the second input argument FIELD is a string specifying an option name, then return the option value VAL corresponding to FIELD from ODE_OPT.

     If called with an optional third input argument, and FIELD is not set in the structure ODE_OPT, then return the default value DEFAULT instead. odeset

 -- STOP_SOLVE =: odeplot (T, Y, FLAG)

     Open a new figure window and plot the solution of an ode problem at each time step during the integration.

     The types and values of the input parameters T and Y depend on the input FLAG that is of type string.  Valid values of FLAG are:

     "init"
          The input T must be a column vector of length 2 with the first and last time step ([TFIRST TLAST].  The input Y contains the initial conditions for the ode problem (Y0).

     ""
          The input T must be a scalar double or vector specifying the time(s) for which the solution in input Y was calculated.

     "done"
          The inputs should be empty, but are ignored if they are present.

     odeplot always returns false, i.e., don't stop the ode solver.

     Example: solve an anonymous implementation of the "Van der Pol" equation and display the results while solving.

          fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];

          opt = odeset ("OutputFcn", @odeplot, "RelTol", 1e-6);
          sol = ode45 (fvdp, [0 20], [2 0], opt);

     Background Information: This function is called by an ode solver function if it was specified in the "OutputFcn" property of an options structure created with odeset.  The ode solver will initially call the function with the syntax odeplot ([TFIRST, TLAST], Y0, "init").  The function initializes internal variables, creates a new figure window, and sets the x limits of the plot.  Subsequently, at each time step during the integration the ode solver calls odeplot (T, Y, []). At the end of the solution the ode solver calls odeplot ([], [], "done") so that odeplot can perform any clean-up actions required. odeset, odeget, ode23, ode45