## Copyright (C) 2014, Jacopo Corno <jacopo.corno@gmail.com>
## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it>
## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net>
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init})
## @deftypefnx {Function File} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init}, @var{opt})
## @deftypefnx {Function File} {[@var{t}, @var{y}] =} ode45 (@dots{}, @var{par1}, @var{par2}, @dots{})
## @deftypefnx {Function File} {[@var{t}, @var{y}, @var{xe}, @var{ye}, @var{ie}] =} ode45 (@dots{})
## @deftypefnx {Function File} {@var{sol} =} ode45 (@var{fun}, @var{trange}, @var{init}, @dots{})
##
## Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs)
## with the well known explicit Dormand-Prince method of order 4.
##
## The first input argument must be a function handle or inline function that
## defines the ODE: @code{y' = f(t,y)}.  The function must accept two inputs
## where the first is time @var{t} and the second is a column vector of
## unknowns @var{y}.
##
## @var{trange} specifies the time interval over which the ODE will be
## evaluated.  Usually, it is a two-element vector specifying the initial and
## final times (@code{[tinit, tfinal]}).  If there are more than two elements
## then the solution will also be evaluated at these intermediate time
## instances unless the integrate function called is
## @command{integrate_n_steps}.  If there is only one time value, then
## @code{ode45} will raise an error unless the options structure has
## non-empty fields named @var{"TimeStepNumber"} and @var{"TimeStepSize"}. 
## If the option @var{"TimeStepSize"} is not empty, then the stepper called
## will be @command{integrate_const}.  If @var{"TimeStepNumber"} is also
## specified then the integrate function @command{integrate_n_steps} will be
## used; otherwise, @command{integrate_adaptive} is used.  For this last
## possibility the user can set the tolerance for the timestep computation by
## changing the option @var{"Tau"}, that has a default value of @math{1e-6}.
##
## The third input argument @var{init} contains the initial value for the
## unknowns.  If this is a row vector then the solution @var{y} will be a matrix
## in which each column is the solution for the corresponding initial value
## in @var{init}.
##
## If present, the fourth input argument specifies options to the ODE solver.
## It is a structure typically generated by @code{odeset}.
##
## The function usually produces just two outputs.  Variable @var{t} is a
## column vector and contains the times where the solution was found.  The
## output @var{y} is a matrix in which each column refers to a different
## unknown of the problem and each row corresponds to a time in @var{t}.
##
## For example, solve an anonymous implementation of the Van der Pol equation
##
## @example
## @group
## fvdp = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
## [T,Y] = ode45 (fvdp, [0 20], [2 0]);
## @end group
## @end example
## @seealso{odeset, odeget}
## @end deftypefn

function varargout = ode45 (vfun, vtrange, vinit, varargin)

  vorder = 5;  # runge_kutta_45_dorpri uses local extrapolation
  vsolver = "ode45";

  ## FIXME: Octave core function's usually do print_usage ()
  ##        rather than displaying full help string when improperly called.
  if (nargin == 0)  # Check number and types of all input arguments
    help (vsolver);
    error ("OdePkg:InvalidArgument",
           "number of input arguments must be greater than zero");
  endif

  if (nargin < 3)
    print_usage ();
  endif
  
  if (nargin >= 4)
    if (! isstruct (varargin{1}))
      ## varargin{1:len} are parameters for vfun
      vodeoptions = odeset;
      vodeoptions.vfunarguments = varargin;
    elseif (length (varargin) > 1)
      ## varargin{1} is an OdePkg options structure vopt
      vodeoptions = odepkg_structure_check (varargin{1}, "ode45");
      vodeoptions.vfunarguments = {varargin{2:length(varargin)}};
    else  # if (isstruct (varargin{1}))
      vodeoptions = odepkg_structure_check (varargin{1}, "ode45");
      vodeoptions.vfunarguments = {};
    endif
  else  # nargin == 3
    vodeoptions = odeset; 
    vodeoptions.vfunarguments = {};
  endif

  if (! isvector (vtrange) || ! isnumeric (vtrange))
    error ("OdePkg:InvalidArgument",
           "second input argument must be a valid vector");
  endif

  if (length (vtrange) < 2
      && (isempty (vodeoptions.TimeStepSize)
          || isempty (vodeoptions.TimeStepNumber)))
    error ("OdePkg:InvalidArgument",
           "second input argument must be a valid vector");
  elseif (vtrange(2) == vtrange(1))
    error ("OdePkg:InvalidArgument",
           "second input argument must be a valid vector");
  else
    vodeoptions.vdirection = sign (vtrange(2) - vtrange(1));
  endif
  vtrange = vtrange(:);

  if (! isvector (vinit) || ! isnumeric (vinit))
    error ("OdePkg:InvalidArgument",
           "third input argument must be a valid numerical value");
  endif
  vinit = vinit(:);

  if (! (isa (vfun, "function_handle")))
    error ("OdePkg:InvalidArgument",
           "first input argument must be a valid function handle");
  endif

  ## Start preprocessing, have a look which options are set in vodeoptions,
  ## check if an invalid or unused option is set
  if (isempty (vodeoptions.TimeStepNumber)
      && isempty (vodeoptions.TimeStepSize))
    integrate_func = "adaptive";
    vodeoptions.vstepsizefixed = false;
  elseif (! isempty (vodeoptions.TimeStepNumber)
          && ! isempty (vodeoptions.TimeStepSize))
    integrate_func = "n_steps";
    vodeoptions.vstepsizefixed = true;
    if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection)
      warning ("OdePkg:InvalidArgument",
               "option 'TimeStepSize' has a wrong sign",
               "it will be corrected automatically");
      vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
    endif
  elseif (isempty (vodeoptions.TimeStepNumber)
          && ! isempty (vodeoptions.TimeStepSize))
    integrate_func = "const";
    vodeoptions.vstepsizefixed = true;
    if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection)
      warning ("OdePkg:InvalidArgument",
               "option 'TimeStepSize' has a wrong sign",
               "it will be corrected automatically");
      vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
    endif
  else
    warning ("OdePkg:InvalidArgument",
             "assuming an adaptive integrate function");
    integrate_func = "adaptive";
  endif

  ## Get the default options that can be set with "odeset" temporarily
  vodetemp = odeset;

  ## Implementation of the option RelTol has been finished.
  ## This option can be set by the user to another value than default value.
  if (isempty (vodeoptions.RelTol) && ! vodeoptions.vstepsizefixed)
    vodeoptions.RelTol = 1e-3;
    warning ("OdePkg:InvalidArgument",
             "Option 'RelTol' not set, new value %f is used",
             vodeoptions.RelTol);
  elseif (! isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed)
    warning ("OdePkg:InvalidArgument",
             "Option 'RelTol' will be ignored if fixed time stamps are given");
  endif

  ## Implementation of the option AbsTol has been finished.
  ## This option can be set by the user to another value than default value.
  if (isempty (vodeoptions.AbsTol) && ! vodeoptions.vstepsizefixed)
    vodeoptions.AbsTol = 1e-6;
    warning ("OdePkg:InvalidArgument",
             "Option 'AbsTol' not set, new value %f is used",
             vodeoptions.AbsTol);
  elseif (! isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed)
    warning ("OdePkg:InvalidArgument",
             "Option 'AbsTol' will be ignored if fixed time stamps are given");
  else
    vodeoptions.AbsTol = vodeoptions.AbsTol(:); # Create column vector
  endif

  ## Implementation of the option NormControl has been finished.
  ## This option can be set by the user to another value than default value.
  if (strcmp (vodeoptions.NormControl, "on"))
    vodeoptions.vnormcontrol = true;
  else 
    vodeoptions.vnormcontrol = false; 
  endif

  ## Implementation of the option NonNegative has been finished.
  ## This option can be set by the user to another value than default value.
  if (! isempty (vodeoptions.NonNegative))
    if (isempty (vodeoptions.Mass))
      vodeoptions.vhavenonnegative = true;
    else
      vodeoptions.vhavenonnegative = false;
      warning ("OdePkg:InvalidArgument",
               "Option 'NonNegative' will be ignored if mass matrix is set");
    endif
  else 
    vodeoptions.vhavenonnegative = false;
  endif

  ## Implementation of the option OutputFcn has been finished.
  ## This option can be set by the user to another value than default value.
  if (isempty (vodeoptions.OutputFcn) && nargout == 0)
    vodeoptions.OutputFcn = @odeplot;
    vodeoptions.vhaveoutputfunction = true;
  elseif (isempty (vodeoptions.OutputFcn))
    vodeoptions.vhaveoutputfunction = false;
  else 
    vodeoptions.vhaveoutputfunction = true;
  endif

  ## Implementation of the option OutputSel has been finished.
  ## This option can be set by the user to another value than default value.
  if (! isempty (vodeoptions.OutputSel))
    vodeoptions.vhaveoutputselection = true;
  else 
    vodeoptions.vhaveoutputselection = false; 
  endif

  ## Implementation of the option OutputSave has been finished.
  ## This option can be set by the user to another value than default value.
  if (isempty (vodeoptions.OutputSave))
    vodeoptions.OutputSave = 1;
  endif

  ## Implementation of the option Refine has been finished.
  ## This option can be set by the user to another value than default value.
  if (vodeoptions.Refine > 0)
    vodeoptions.vhaverefine = true;
  else 
    vodeoptions.vhaverefine = false;
  endif

  ## Implementation of the option Stats has been finished.
  ## This option can be set by the user to another value than default value.

  ## Implementation of the option InitialStep has been finished.
  ## This option can be set by the user to another value than default value.
  if (isempty (vodeoptions.InitialStep) && strcmp (integrate_func, "adaptive"))
    vodeoptions.InitialStep = ...
      vodeoptions.vdirection * starting_stepsize (vorder, vfun, vtrange(1),
                                                  vinit,
                                                  vodeoptions.AbsTol,
                                                  vodeoptions.RelTol,
                                                  vodeoptions.vnormcontrol);
    warning ("OdePkg:InvalidArgument",
             "option 'InitialStep' not set, estimated value %f is used",
             vodeoptions.InitialStep);
  elseif(isempty (vodeoptions.InitialStep))
    vodeoptions.InitialStep = odeget (vodeoptions, "TimeStepSize");
  endif

  ## Implementation of the option MaxStep has been finished. This option
  ## can be set by the user to another value than default value.
  if (isempty (vodeoptions.MaxStep) && ! vodeoptions.vstepsizefixed)
    vodeoptions.MaxStep = abs (vtrange(end) - vtrange(1)) / 10;
    warning ("OdePkg:InvalidArgument",
             "Option 'MaxStep' not set, new value %f is used",
             vodeoptions.MaxStep);
  endif

  ## Implementation of the option Events has been finished.
  ## This option can be set by the user to another value than default value.
  if (! isempty (vodeoptions.Events))
    vodeoptions.vhaveeventfunction = true;
  else 
    vodeoptions.vhaveeventfunction = false;
  endif

  ## The options "Jacobian", "JPattern" and "Vectorized" will be ignored
  ## by this solver because this solver uses an explicit Runge-Kutta method
  ## and therefore no Jacobian calculation is necessary.
  if (! isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
    warning ("OdePkg:InvalidArgument",
             "option 'Jacobian' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.JPattern, vodetemp.JPattern))
    warning ("OdePkg:InvalidArgument",
             "option 'JPattern' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
    warning ("OdePkg:InvalidArgument",
             "option 'Vectorized' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
    warning ("OdePkg:InvalidArgument",
             "option 'NewtonTol' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.MaxNewtonIterations,vodetemp.MaxNewtonIterations))
    warning ("OdePkg:InvalidArgument",
             "option 'MaxNewtonIterations' will be ignored by this solver");
  endif

  ## Implementation of the option Mass has been finished.
  ## This option can be set by the user to another value than default value.
  if (! isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
    vhavemasshandle = false;
    vmass = vodeoptions.Mass;  # constant mass
  elseif (isa (vodeoptions.Mass, "function_handle"))
    vhavemasshandle = true;    # mass defined by a function handle
  else  # no mass matrix - creating a diag-matrix of ones for mass
    vhavemasshandle = false;   # vmass = diag (ones (length (vinit), 1), 0);
  endif

  ## Implementation of the option MStateDependence has been finished.
  ## This option can be set by the user to another value than default value.
  if (strcmp (vodeoptions.MStateDependence, "none"))
    vmassdependence = false;
  else
    vmassdependence = true;
  endif

  ## Other options that are not used by this solver.
  ## Print a warning message to tell the user that the option is ignored.
  if (! isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
    warning ("OdePkg:InvalidArgument",
             "option 'MvPattern' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
    warning ("OdePkg:InvalidArgument",
             "option 'MassSingular' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
    warning ("OdePkg:InvalidArgument",
             "option 'InitialSlope' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
    warning ("OdePkg:InvalidArgument",
             "option 'MaxOrder' will be ignored by this solver");
  endif

  if (! isequal (vodeoptions.BDF, vodetemp.BDF))
    warning ("OdePkg:InvalidArgument",
             "option 'BDF' will be ignored by this solver");
  endif

  ## Starting the initialisation of the core solver ode45
  SubOpts = vodeoptions;
  
  if (vhavemasshandle)   # Handle only the dynamic mass matrix,
    if (vmassdependence) # constant mass matrices have already
      vmass = @(t,x) vodeoptions.Mass (t, x, vodeoptions.vfunarguments{:});
      vfun = @(t,x) vmass (t, x, vodeoptions.vfunarguments{:}) ...
             \ vfun (t, x, vodeoptions.vfunarguments{:});
    else                 # if (vmassdependence == false)
      vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:});
      vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ...
             \ vfun (t, x, vodeoptions.vfunarguments{:});
    endif
  endif

  switch (integrate_func)
    case "adaptive"
      solution = integrate_adaptive (@runge_kutta_45_dorpri,
                                     vorder, vfun, vtrange, vinit, SubOpts);
    case "n_steps"
      solution = integrate_n_steps (@runge_kutta_45_dorpri,
                                    vfun, vtrange(1), vinit,
                                    vodeoptions.TimeStepSize,
                                    vodeoptions.TimeStepNumber, SubOpts);
    case "const"
      solution = integrate_const (@runge_kutta_45_dorpri,
                                  vfun, vtrange, vinit,
                                  vodeoptions.TimeStepSize, SubOpts);
  endswitch

  ## Postprocessing, do whatever when terminating integration algorithm
  if (vodeoptions.vhaveoutputfunction)  # Cleanup plotter
    feval (vodeoptions.OutputFcn, solution.t(end),
           solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
  endif
  if (vodeoptions.vhaveeventfunction)   # Cleanup event function handling
    odepkg_event_handle (vodeoptions.Events, solution.t(end),
                         solution.x(end,:)', "done",
                         vodeoptions.vfunarguments{:});
  endif

  ## Print additional information if option Stats is set
  if (strcmp (vodeoptions.Stats, "on"))
    vhavestats = true;
    vnsteps    = solution.vcntloop-2;        # vcntloop from 2..end
    vnfailed   = (solution.vcntcycles-1)-(vnsteps)+1; # vcntcycl from 1..end
    vnfevals   = 7*(solution.vcntcycles-1);  # number of ode evaluations
    vndecomps  = 0;                          # number of LU decompositions
    vnpds      = 0;                          # number of partial derivatives
    vnlinsols  = 0;                          # no. of linear systems solutions
    ## Print cost statistics if no output argument is given
    if (nargout == 0)
      printf ("Number of successful steps: %d\n", vnsteps);
      printf ("Number of failed attempts:  %d\n", vnfailed);
      printf ("Number of function calls:   %d\n", vnfevals);
    endif
  else
    vhavestats = false;
  endif

  if (nargout == 2)
    varargout{1} = solution.t;      # Time stamps are first output argument
    varargout{2} = solution.x;      # Results are second output argument
  elseif (nargout == 1)
    varargout{1}.x = solution.t;    # Time stamps are saved in field x
    varargout{1}.y = solution.x;    # Results are saved in field y
    varargout{1}.solver = vsolver;  # Solver name is saved in field solver
    if (vodeoptions.vhaveeventfunction) 
      varargout{1}.ie = solution.vevent{2};  # Index info which event occurred
      varargout{1}.xe = solution.vevent{3};  # Time info when an event occurred
      varargout{1}.ye = solution.vevent{4};  # Results when an event occurred
    endif
    if (vhavestats)
      varargout{1}.stats = struct ();
      varargout{1}.stats.nsteps   = vnsteps;
      varargout{1}.stats.nfailed  = vnfailed;
      varargout{1}.stats.nfevals  = vnfevals;
      varargout{1}.stats.npds     = vnpds;
      varargout{1}.stats.ndecomps = vndecomps;
      varargout{1}.stats.nlinsols = vnlinsols;
    endif
  elseif (nargout == 5)
    varargout = cell (1,5);
    varargout{1} = solution.t;
    varargout{2} = solution.x;
    if (vodeoptions.vhaveeventfunction) 
      varargout{3} = solution.vevent{3};     # Time info when an event occurred
      varargout{4} = solution.vevent{4};     # Results when an event occurred
      varargout{5} = solution.vevent{2};     # Index info which event occurred
    endif
  endif

endfunction


## We are using the "Van der Pol" implementation for all tests that are done
## for this function.
## For further tests we also define a reference solution (computed at high
## accuracy)
%!function [ydot] = fpol (vt, vy)  # The Van der Pol
%! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%!endfunction
%!function [vref] = fref ()        # The computed reference solution
%! vref = [0.32331666704577, -1.83297456798624];
%!endfunction
%!function [vjac] = fjac (vt, vy, varargin) # its Jacobian
%! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
%!endfunction
%!function [vjac] = fjcc (vt, vy, varargin) # sparse type
%! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2])
%!endfunction
%!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
%! vval = fpol (vt, vy, varargin); # We use the derivatives
%! vtrm = zeros (2,1);             # that's why component 2
%! vdir = ones (2,1);              # seems to not be exact
%!endfunction
%!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
%! vval = fpol (vt, vy, varargin); # We use the derivatives
%! vtrm = ones (2,1);              # that's why component 2
%! vdir = ones (2,1);              # seems to not be exact
%!endfunction
%!function [vmas] = fmas (vt, vy, varargin)
%! vmas = [1, 0; 0, 1];            # Dummy mass matrix for tests
%!endfunction
%!function [vmas] = fmsa (vt, vy, varargin)
%! vmas = sparse ([1, 0; 0, 1]);   # A sparse dummy matrix
%!endfunction
%!function [vout] = fout (vt, vy, vflag, varargin)
%! if (regexp (char (vflag), 'init') == 1)
%!   if (any (size (vt) != [2, 1])) error ('"fout" step "init"'); endif
%! elseif (isempty (vflag))
%!   if (any (size (vt) != [1, 1])) error ('"fout" step "calc"'); endif
%!   vout = false;
%! elseif (regexp (char (vflag), 'done') == 1)
%!   if (any (size (vt) != [1, 1])) error ('"fout" step "done"'); endif
%! else
%!   error ('"fout" invalid vflag');
%! endif
%!endfunction
%!
%! ## Turn off output of warning messages for all tests, turn them on
%! ## again if the last test is called
%!error  # ouput argument
%! warning ("off", "OdePkg:InvalidArgument");
%! B = ode45 (1, [0 25], [3 15 1]);
%!error  # input argument number one
%! [vt, vy] = ode45 (1, [0 25], [3 15 1]);
%!error  # input argument number two
%! [vt, vy] = ode45 (@fpol, 1, [3 15 1]);
%!test  # two output arguments
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
%!test  # not too many steps
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
%! assert (size (vt) < 20);
%!test  # anonymous function instead of real function
%! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%! [vt, vy] = ode45 (fvdb, [0 2], [2 0]);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
%!test  # extra input arguments passed through
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], 12, 13, "KL");
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
%!test  # empty OdePkg structure *but* extra input arguments
%! vopt = odeset;
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt, 12, 13, "KL");
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
%!error  # strange OdePkg structure
%! vopt = struct ("foo", 1);
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
%!test  # Solve vdp in fixed step sizes
%! vopt = odeset("TimeStepSize", 0.1);
%! [vt, vy] = ode45 (@fpol, [0,2], [2 0], vopt);
%! assert (vt(:), [0:0.1:2]', 1e-2);
%!test  # Solve another anonymous function below zero
%! vref = [0, 14.77810590694212];
%! [vt, vy] = ode45 (@(t,y) y, [-2 0], 2);
%! assert ([vt(end), vy(end,:)], vref, 1e-1);
%!test  # InitialStep option
%! vopt = odeset ("InitialStep", 1e-8);
%! [vt, vy] = ode45 (@fpol, [0 0.2], [2 0], vopt);
%! assert ([vt(2)-vt(1)], [1e-8], 1e-9);
%!test  # MaxStep option
%! vopt = odeset ("MaxStep", 1e-3);
%! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt);
%! assert ([vsol.x(5)-vsol.x(4)], [1e-3], 1e-3);
%!test  # Solve with intermidiate step
%! vsol = ode45 (@fpol, [0 1 2], [2 0]);
%! assert (any((vsol.x-1) == 0));
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # Solve in backward direction starting at t=0
%! vref = [-1.205364552835178, 0.951542399860817];
%! vsol = ode45 (@fpol, [0 -2], [2 0]);
%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
%!test  # Solve in backward direction starting at t=2
%! vref = [-1.205364552835178, 0.951542399860817];
%! vsol = ode45 (@fpol, [2 -2], fref);
%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
%!test  # Solve in backward direction starting at t=2, with intermidiate step
%! vref = [-1.205364552835178, 0.951542399860817];
%! vsol = ode45 (@fpol, [2 0 -2], fref);
%! idx = find(vsol.x < 0, 1, "first") - 1;
%! assert ([vsol.x(idx), vsol.y(idx,:)], [0 2 0], 1e-2);
%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
%!test  # Solve another anonymous function in backward direction
%! vref = [-1, 0.367879437558975];
%! vsol = ode45 (@(t,y) y, [0 -1], 1);
%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
%!test  # Solve another anonymous function below zero
%! vref = [0, 14.77810590694212];
%! vsol = ode45 (@(t,y) y, [-2 0], 2);
%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
%!test  # Solve in backward direction starting at t=0 with MaxStep option
%! vref = [-1.205364552835178, 0.951542399860817];
%! vopt = odeset ("MaxStep", 1e-3);
%! vsol = ode45 (@fpol, [0 -2], [2 0], vopt);
%! assert ([abs(vsol.x(8)-vsol.x(7))], [1e-3], 1e-3);
%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
%!test  # AbsTol option
%! vopt = odeset ("AbsTol", 1e-5);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # AbsTol and RelTol option
%! vopt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # RelTol and NormControl option -- higher accuracy
%! vopt = odeset ("RelTol", 1e-8, "NormControl", "on");
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-5);
%!test  # Keeps initial values while integrating
%! vopt = odeset ("NonNegative", 2);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5);
%!test  # Details of OutputSel and Refine can't be tested
%! vopt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%!test  # Details of OutputSave can't be tested
%! vopt = odeset ("OutputSave", 1, "OutputSel", 1);
%! vsla = ode45 (@fpol, [0 2], [2 0], vopt);
%! vopt = odeset ("OutputSave", 2);
%! vslb = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert (length (vsla.x) + 1 >= 2 * length (vslb.x))
%!test  # Stats must add further elements in vsol
%! vopt = odeset ("Stats", "on");
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert (isfield (vsol, "stats"));
%! assert (isfield (vsol.stats, "nsteps"));
%!test  # Events option add further elements in vsol
%! vopt = odeset ("Events", @feve);
%! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
%! assert (isfield (vsol, "ie"));
%! assert (vsol.ie(1), 2);
%! assert (isfield (vsol, "xe"));
%! assert (isfield (vsol, "ye"));
%!test  # Events option, now stop integration
%! vopt = odeset ("Events", @fevn, "NormControl", "on");
%! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
%! assert ([vsol.ie, vsol.xe, vsol.ye], 
%!         [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
%!test  # Events option, five output arguments
%! vopt = odeset ("Events", @fevn, "NormControl", "on");
%! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 10], [2 0], vopt);
%! assert ([vie, vxe, vye],
%!         [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
%!test  # Jacobian option
%! vopt = odeset ("Jacobian", @fjac);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # Jacobian option and sparse return value
%! vopt = odeset ("Jacobian", @fjcc);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!
%!## test for JPattern option is missing
%!## test for Vectorized option is missing
%!## test for NewtonTol option is missing
%!## test for MaxNewtonIterations option is missing
%!
%!test  # Mass option as function
%! vopt = odeset ("Mass", @fmas);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # Mass option as matrix
%! vopt = odeset ("Mass", eye (2,2));
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # Mass option as sparse matrix
%! vopt = odeset ("Mass", sparse (eye (2,2)));
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # Mass option as function and sparse matrix
%! vopt = odeset ("Mass", @fmsa);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # Mass option as function and MStateDependence
%! vopt = odeset ("Mass", @fmas, "MStateDependence", "strong");
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test  # Set BDF option to something else than default
%! vopt = odeset ("BDF", "on");
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
%!
%!## test for MvPattern option is missing
%!## test for InitialSlope option is missing
%!## test for MaxOrder option is missing
%!
%! warning ("on", "OdePkg:InvalidArgument");