— Command: [] = odebdi (@fun, slot, y0, dy0, [opt], [P1, P2, ...])
— Command: [sol] = odebdi (@fun, slot, y0, dy0, [opt], [P1, P2, ...])
— Command: [t, y, [xe, ye, ie]] = odebdi (@fun, slot, y0, dy0, [opt], [P1, P2, ...])

This function file can be used to solve a set of stiff implicit differential equations (IDEs). This function file is a wrapper file that uses Jeff Cash's Fortran solver mebdfi.f.

If this function is called with no return argument then plot the solution over time in a figure window while solving the set of IDEs that are defined in a function and specified by the function handle @fun. The second input argument slot is a double vector that defines the time slot, y0 is a double vector that defines the initial values of the states, dy0 is a double vector that defines the initial values of the derivatives, opt can optionally be a structure array that keeps the options created with the command odeset and par1, par2, ... can optionally be other input arguments of any type that have to be passed to the function defined by @fun.

If this function is called with one return argument then return the solution sol of type structure array after solving the set of IDEs. The solution sol has the fields x of type double column vector for the steps chosen by the solver, y of type double column vector for the solutions at each time step of x, solver of type string for the solver name and optionally the extended time stamp information xe, the extended solution information ye and the extended index information ie all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument opt.

If this function is called with more than one return argument then return the time stamps t, the solution values y and optionally the extended time stamp information xe, the extended solution information ye and the extended index information ie all of type double column vector.

For example,

          function res = odepkg_equations_ilorenz (t, y, yd)
            res = [10 * (y(2) - y(1)) - yd(1);
                   y(1) * (28 - y(3)) - yd(2);
                   y(1) * y(2) - 8/3 * y(3) - yd(3)];
          endfunction
          
          vopt = odeset ("InitialStep\", 1e-3, \"MaxStep\", 1e-1, \\
                         \"OutputFcn\", @odephas3, \"Refine\", 5);
          odebdi (@odepkg_equations_ilorenz, [0, 25], [3 15 1], \\
                  [120 81 42.333333], vopt);