# HG changeset patch
# User Stefan Miereis <stefan.miereis@gmx.de>
# Date 1449941606 -3600
#      Sat Dec 12 18:33:26 2015 +0100
# Node ID b0c8bc80acd46a3212b2290e04d5764782bf0d0a
# Parent  61892a245653c0a09f5cc4697e591c9d81d2f089
Code Sprint 2015: Move ode23 from odepkg to core

diff -r 61892a245653 -r b0c8bc80acd4 doc/interpreter/diffeq.txi
--- a/doc/interpreter/diffeq.txi	Sat Dec 12 08:31:22 2015 -0800
+++ b/doc/interpreter/diffeq.txi	Sat Dec 12 18:33:26 2015 +0100
@@ -143,12 +143,17 @@
     method.  This is a fourth--order accurate integrator therefore the local
     error normally expected is @math{O(h^5)}.  This solver requires six
     function evaluations per integration step.
+
+    @item @code{ode23} Integrates a system of non--stiff ordinary differential equations
+    (non--stiff ODEs and DAEs) with the well known explicit Runge--Kutta method of order (2,3)
   @end itemize
 @end itemize
 
 
 @DOCSTRING(ode45)
 
+@DOCSTRING(ode23)
+
 @DOCSTRING(odeset)
 
 @DOCSTRING(odeget)
diff -r 61892a245653 -r b0c8bc80acd4 scripts/ode/module.mk
--- a/scripts/ode/module.mk	Sat Dec 12 08:31:22 2015 -0800
+++ b/scripts/ode/module.mk	Sat Dec 12 18:33:26 2015 +0100
@@ -17,6 +17,7 @@
 
 scripts_ode_FCN_FILES = \
   scripts/ode/ode45.m \
+  scripts/ode/ode23.m \
   scripts/ode/odeset.m \
   scripts/ode/odeget.m
 
diff -r 61892a245653 -r b0c8bc80acd4 scripts/ode/ode23.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/ode/ode23.m	Sat Dec 12 18:33:26 2015 +0100
@@ -0,0 +1,747 @@
+## Copyright (C) 2006-2015 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/>.
+
+## Author: Thomas Treichl <treichl@users.sourceforge.net>
+
+%# -*- texinfo -*-
+%# @deftypefn  {Function File} {[@var{}] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
+%# @deftypefnx {Command} {[@var{sol}] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
+%# @deftypefnx {Command} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
+%#
+%# This function file can be used to solve a set of non--stiff ordinary differential equations (non--stiff ODEs) or non--stiff differential algebraic equations (non--stiff DAEs) with the well known explicit Runge--Kutta method of order (2,3).
+%#
+%# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of ODEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
+%#
+%# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of ODEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{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 @var{opt}.
+%#
+%# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
+%#
+%# For example, solve an anonymous implementation of the Van der Pol equation
+%#
+%# @example
+%# fvdb = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
+%#
+%# opt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, \
+%#          "NormControl", "on", "OutputFcn", @@odeplot);
+%# ode23 (fvdb, [0 20], [2 0], opt);
+%# @end example
+%# @end deftypefn
+%#
+%# @seealso{ode45, odeset, odeget}
+
+function [varargout] = ode23 (fun, slot, init, varargin)
+
+  if (nargin == 0) # Check number and types of all input arguments
+    help ("ode23");
+    error ("Octave:InvalidArgument", ...
+      "ode23: Number of input arguments must be greater than zero");
+
+  elseif (nargin < 3)
+    print_usage;
+
+  elseif !(isa (fun, "function_handle") || isa (fun, "inline"))
+    error ("Octave:InvalidArgument", ...
+      "ode23: First input argument must be a valid function handle");
+
+  elseif (!isvector (slot) || length (slot) < 2)
+    error ("Octave:InvalidArgument", ...
+      "ode23: Second input argument must be a valid vector");
+
+  elseif (!isvector (init) || !isnumeric (init))
+    error ("Octave:InvalidArgument", ...
+      "ode23: Third input argument must be a valid numerical value");
+
+  elseif (nargin >= 4)
+
+    if (!isstruct (varargin{1}))
+      ## varargin{1:len} are parameters for fun
+      odeoptions = odeset;
+      funarguments = varargin;
+
+    elseif (length (varargin) > 1)
+      ## varargin{1} is an OdePkg options structure opt
+      odeoptions = odepkg_structure_check (varargin{1}, "ode23");
+      funarguments = {varargin{2:length(varargin)}};
+
+    else # if (isstruct (varargin{1}))
+      odeoptions = odepkg_structure_check (varargin{1}, "ode23");
+      funarguments = {};
+
+    endif
+
+  else # if (nargin == 3)
+    odeoptions = odeset;
+    funarguments = {};
+  endif
+
+  ## Start preprocessing, have a look which options are set in
+  ## odeoptions, check if an invalid or unused option is set
+  slot = slot(:).';     # Create a row vector
+  init = init(:).';     # Create a row vector
+  if (length (slot) > 2) # Step size checking
+    stepsizefixed = true;
+  else
+    stepsizefixed = false;
+  endif
+
+  ## Get the default options that can be set with "odeset" temporarily
+  odetemp = 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 (odeoptions.RelTol) && !stepsizefixed)
+    odeoptions.RelTol = 1e-6;
+    warning ("Octave:InvalidArgument", ...
+      "Option \"RelTol\" not set, new value %f is used", odeoptions.RelTol);
+  elseif (!isempty (odeoptions.RelTol) && stepsizefixed)
+    warning ("Octave: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 (odeoptions.AbsTol) && !stepsizefixed)
+    odeoptions.AbsTol = 1e-6;
+    warning ("Octave:InvalidArgument", ...
+      "Option \"AbsTol\" not set, new value %f is used", odeoptions.AbsTol);
+  elseif (!isempty (odeoptions.AbsTol) && stepsizefixed)
+    warning ("Octave:InvalidArgument", ...
+      "Option \"AbsTol\" will be ignored if fixed time stamps are given");
+  else
+    odeoptions.AbsTol = odeoptions.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 (odeoptions.NormControl, "on")) normcontrol = true;
+  else normcontrol = 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 (odeoptions.NonNegative))
+    if (isempty (odeoptions.Mass)), havenonnegative = true;
+    else
+      havenonnegative = false;
+      warning ("Octave:InvalidArgument", ...
+        "Option \"NonNegative\" will be ignored if mass matrix is set");
+    endif
+  else havenonnegative = 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 (odeoptions.OutputFcn) && nargout == 0)
+    odeoptions.OutputFcn = @odeplot;
+    haveoutputfunction = true;
+  elseif (isempty (odeoptions.OutputFcn)), haveoutputfunction = false;
+  else haveoutputfunction = 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 (odeoptions.OutputSel)), haveoutputselection = true;
+  else haveoutputselection = 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 (odeoptions.OutputSave)), odeoptions.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 (odeoptions.Refine > 0), haverefine = true;
+  else haverefine = 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 (odeoptions.InitialStep) && !stepsizefixed)
+    odeoptions.InitialStep = (slot(1,2) - slot(1,1)) / 10;
+    odeoptions.InitialStep = odeoptions.InitialStep / 10^odeoptions.Refine;
+    warning ("Octave:InvalidArgument", ...
+      "Option \"InitialStep\" not set, new value %f is used", odeoptions.InitialStep);
+  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 (odeoptions.MaxStep) && !stepsizefixed)
+    odeoptions.MaxStep = abs (slot(1,2) - slot(1,1)) / 10;
+    warning ("Octave:InvalidArgument", ...
+      "Option \"MaxStep\" not set, new value %f is used", odeoptions.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 (odeoptions.Events)), haveeventfunction = true;
+  else haveeventfunction = 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 (odeoptions.Jacobian, odetemp.Jacobian))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"Jacobian\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.JPattern, odetemp.JPattern))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"JPattern\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.Vectorized, odetemp.Vectorized))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"Vectorized\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.NewtonTol, odetemp.NewtonTol))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"NewtonTol\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.MaxNewtonIterations,...
+                odetemp.MaxNewtonIterations))
+    warning ("Octave: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 (odeoptions.Mass) && isnumeric (odeoptions.Mass))
+    havemasshandle = false; mass = odeoptions.Mass; %# constant mass
+  elseif (isa (odeoptions.Mass, "function_handle"))
+    havemasshandle = true; # mass defined by a function handle
+  else # no mass matrix - creating a diag-matrix of ones for mass
+    havemasshandle = false; # mass = diag (ones (length (init), 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 (odeoptions.MStateDependence, "none"))
+    massdependence = false;
+  else massdependence = true;
+  endif
+
+  ## Other options that are not used by this solver. Print a warning
+  ## message to tell the user that the option(s) is/are ignored.
+  if (!isequal (odeoptions.MvPattern, odetemp.MvPattern))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"MvPattern\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.MassSingular, odetemp.MassSingular))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"MassSingular\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.InitialSlope, odetemp.InitialSlope))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"InitialSlope\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.MaxOrder, odetemp.MaxOrder))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"MaxOrder\" will be ignored by this solver");
+  endif
+  if (!isequal (odeoptions.BDF, odetemp.BDF))
+    warning ("Octave:InvalidArgument", ...
+      "Option \"BDF\" will be ignored by this solver");
+  endif
+
+  ## Starting the initialisation of the core solver ode23
+  timestamp  = slot(1,1);           # timestamp = start time
+  timelength = length (slot);       # length needed if fixed steps
+  timestop   = slot(1,timelength); # stop time = last value
+  ## 20110611, reported by Nils Strunk
+  ## Make it possible to solve equations from negativ to zero,
+  ## eg. res = ode23 (@(t,y) y, [-2 0], 2);
+  direction  = sign (timestop - timestamp); %# Direction flag
+
+  if (!stepsizefixed)
+    if (sign (odeoptions.InitialStep) == direction)
+      stepsize = odeoptions.InitialStep;
+    else # Fix wrong direction of InitialStep.
+      stepsize = - odeoptions.InitialStep;
+    endif
+    minstepsize = (timestop - timestamp) / (1/eps);
+  else # If step size is given then use the fixed time steps
+    stepsize = slot(1,2) - slot(1,1);
+    minstepsize = sign (stepsize) * eps;
+  endif
+
+  retvaltime = timestamp; # first timestamp output
+  retvalresult = init;    # first solution output
+
+  ## Initialize the OutputFcn
+  if (haveoutputfunction)
+    if (haveoutputselection) retout = retvalresult(odeoptions.OutputSel);
+    else retout = retvalresult; endif
+    feval (odeoptions.OutputFcn, slot.', ...
+      retout.', "init", funarguments{:});
+  endif
+
+  ## Initialize the EventFcn
+  if (haveeventfunction)
+    odepkg_event_handle (odeoptions.Events, timestamp, ...
+      retvalresult.', "init", funarguments{:});
+  endif
+
+  pow = 1/3;            # 20071016, reported by Luis Randez
+  a = [  0, 0, 0;       # The Runge-Kutta-Fehlberg 2(3) coefficients
+        1/2, 0, 0;       # Coefficients proved on 20060827
+         -1, 2, 0];      # See p.91 in Ascher & Petzold
+  b2 = [0; 1; 0];       # 2nd and 3rd order
+  b3 = [1/6; 2/3; 1/6]; # b-coefficients
+  c = sum (a, 2);
+
+  ## The solver main loop - stop if the endpoint has been reached
+  cntloop = 2; cntcycles = 1; vu = init; vk = vu.' * zeros(1,3);
+  cntiter = 0; unhandledtermination = true; cntsave = 2;
+  while ((direction * (timestamp) < direction * (timestop)) && ...
+         (direction * (stepsize) >= direction * (minstepsize)))
+
+    ## Hit the endpoint of the time slot exactely
+    if (direction * (timestamp + stepsize) > direction * timestop)
+      ## stepsize = timestop - direction * timestamp;
+      ## 20110611, reported by Nils Strunk
+      ## The endpoint of the time slot must be hit exactly,
+      ## eg. sol = ode23 (@(t,y) y, [0 -1], 1);
+      stepsize = direction * abs (abs (timestop) - abs (timestamp));
+    endif
+
+    ## Estimate the three results when using this solver
+    for j = 1:3
+      thetime  = timestamp + c(j,1) * stepsize;
+      theinput = vu.' + stepsize * vk(:,1:j-1) * a(j,1:j-1).';
+      if (havemasshandle)   # Handle only the dynamic mass matrix,
+        if (massdependence) # constant mass matrices have already
+          mass = feval ...  # been set before (if any)
+            (odeoptions.Mass, thetime, theinput, funarguments{:});
+        else                 # if (massdependence == false)
+          mass = feval ...  # then we only have the time argument
+            (odeoptions.Mass, thetime, funarguments{:});
+        endif
+        vk(:,j) = mass \ feval ...
+          (fun, thetime, theinput, funarguments{:});
+      else
+        vk(:,j) = feval ...
+          (fun, thetime, theinput, funarguments{:});
+      endif
+    endfor
+
+    ## Compute the 2nd and the 3rd order estimation
+    y2 = vu.' + stepsize * (vk * b2);
+    y3 = vu.' + stepsize * (vk * b3);
+    if (havenonnegative)
+      vu(odeoptions.NonNegative) = abs (vu(odeoptions.NonNegative));
+      y2(odeoptions.NonNegative) = abs (y2(odeoptions.NonNegative));
+      y3(odeoptions.NonNegative) = abs (y3(odeoptions.NonNegative));
+    endif
+    if (haveoutputfunction && haverefine)
+      saveVUForRefine = vu;
+    endif
+
+    ## Calculate the absolute local truncation error and the acceptable error
+    if (!stepsizefixed)
+      if (!normcontrol)
+        delta = abs (y3 - y2);
+        tau = max (odeoptions.RelTol * abs (vu.'), odeoptions.AbsTol);
+      else
+        delta = norm (y3 - y2, Inf);
+        tau = max (odeoptions.RelTol * max (norm (vu.', Inf), 1.0), ...
+                    odeoptions.AbsTol);
+      endif
+    else # if (stepsizefixed == true)
+      delta = 1; tau = 2;
+    endif
+
+    ## If the error is acceptable then update the retval variables
+    if (all (delta <= tau))
+      timestamp = timestamp + stepsize;
+      vu = y3.'; # MC2001: the higher order estimation as "local extrapolation"
+      ## Save the solution every odeoptions.OutputSave steps
+      if (mod (cntloop-1,odeoptions.OutputSave) == 0)
+        retvaltime(cntsave,:) = timestamp;
+        retvalresult(cntsave,:) = vu;
+        cntsave = cntsave + 1;
+      endif
+      cntloop = cntloop + 1; cntiter = 0;
+
+      ## Call plot only if a valid result has been found, therefore this
+      ## code fragment has moved here. Stop integration if plot function
+      ## returns false
+      if (haveoutputfunction)
+        for cnt = 0:odeoptions.Refine # Approximation between told and t
+          if (haverefine)              # Do interpolation
+            approxtime = (cnt + 1) * stepsize / (odeoptions.Refine + 2);
+            approxvals = saveVUForRefine.' + approxtime * (vk * b3);
+            approxtime = (timestamp - stepsize) + approxtime;
+          else
+            approxvals = vu.';
+            approxtime = timestamp;
+          endif
+          if (haveoutputselection)
+            approxvals = approxvals(odeoptions.OutputSel);
+          endif
+          pltret = feval (odeoptions.OutputFcn, approxtime, ...
+            approxvals, [], funarguments{:});
+          if pltret # Leave refinement loop
+            break;
+          endif
+        endfor
+        if (pltret) # Leave main loop
+          unhandledtermination = false;
+          break;
+        endif
+      endif
+
+      ## Call event only if a valid result has been found, therefore this
+      ## code fragment has moved here. Stop integration if veventbreak is
+      ## true
+      if (haveeventfunction)
+        event = ...
+          odepkg_event_handle (odeoptions.Events, timestamp, ...
+            vu(:), [], funarguments{:});
+        if (!isempty (event{1}) && event{1} == 1)
+          retvaltime(cntloop-1,:) = event{3}(end,:);
+          retvalresult(cntloop-1,:) = event{4}(end,:);
+          unhandledtermination = false; break;
+        endif
+      endif
+    endif # If the error is acceptable ...
+
+    ## Update the step size for the next integration step
+    if (!stepsizefixed)
+      ## 20080425, reported by Marco Caliari
+      ## delta cannot be negative (because of the absolute value that
+      ## has been introduced) but it could be 0, then replace the zeros
+      ## with the maximum value of delta
+      delta(find (delta == 0)) = max (delta);
+      ## It could happen that max (delta) == 0 (ie. that the original
+      ## delta was 0), in that case we double the previous stepsize
+      delta(find (delta == 0)) = max (tau) .* (0.4 ^ (1 / pow));
+
+      if (direction == 1)
+        stepsize = min (odeoptions.MaxStep, ...
+           min (0.8 * stepsize * (tau ./ delta) .^ pow));
+      else
+        stepsize = max (- odeoptions.MaxStep, ...
+          max (0.8 * stepsize * (tau ./ delta) .^ pow));
+      endif
+
+    else # if (stepsizefixed)
+      if (cntloop <= timelength)
+        stepsize = slot(cntloop) - slot(cntloop-1);
+      else # Get out of the main integration loop
+        break;
+      endif
+    endif
+
+    ## Update counters that count the number of iteration cycles
+    cntcycles = cntcycles + 1; # Needed for cost statistics
+    cntiter = cntiter + 1;     # Needed to find iteration problems
+
+    ## Stop solving because the last 1000 steps no successful valid
+    ## value has been found
+    if (cntiter >= 5000)
+      error (["ode23: Solving has not been successful. The iterative", ...
+        " integration loop exited at time t = %f before endpoint at", ...
+        " tend = %f was reached. This happened because the iterative", ...
+        " integration loop does not find a valid solution at this time", ...
+        " stamp. Try to reduce the value of \"InitialStep\" and/or", ...
+        " \"MaxStep\" with the command \"odeset\".\n"], timestamp, timestop);
+    endif
+
+  endwhile # The main loop
+
+  ## Check if integration of the ode has been successful
+  if (direction * timestamp < direction * timestop)
+    if (unhandledtermination == true)
+      error ("Octave:InvalidArgument", ...
+        ["ode23: Solving has not been successful. The iterative", ...
+         " integration loop exited at time t = %f", ...
+         " before endpoint at tend = %f was reached. This may", ...
+         " happen if the stepsize grows smaller than defined in", ...
+         " minstepsize. Try to reduce the value of \"InitialStep\" and/or", ...
+         " \"MaxStep\" with the command \"odeset\".\n"], timestamp, timestop);
+    else
+      warning ("Octave:InvalidArgument", ...
+        ["ode23: Solver has been stopped by a call of \"break\" in", ...
+         " the main iteration loop at time t = %f before endpoint at", ...
+         " tend = %f was reached. This may happen because the @odeplot", ...
+         " function returned \"true\" or the @event function returned \"true\"."], ...
+         timestamp, timestop);
+    endif
+  endif
+
+  ## Postprocessing, do whatever when terminating integration algorithm
+  if (haveoutputfunction) # Cleanup plotter
+    feval (odeoptions.OutputFcn, timestamp, ...
+      vu.', "done", funarguments{:});
+  endif
+  if (haveeventfunction)  # Cleanup event function handling
+    odepkg_event_handle (odeoptions.Events, timestamp, ...
+      vu.', "done", funarguments{:});
+  endif
+  ## Save the last step, if not already saved
+  if (mod (cntloop-2,odeoptions.OutputSave) != 0)
+    retvaltime(cntsave,:) = timestamp;
+    retvalresult(cntsave,:) = vu;
+  endif
+
+  ## Print additional information if option Stats is set
+  if (strcmp (odeoptions.Stats, 'on'))
+    havestats = true;
+    nsteps    = cntloop-2;                     # cntloop from 2..end
+    nfailed   = (cntcycles-1)-(cntloop-2)+1;   # cntcycl from 1..end
+    nfevals   = 3*(cntcycles-1);               # number of ode evaluations
+    ndecomps  = 0;                             # number of LU decompositions
+    npds      = 0;                             # number of partial derivatives
+    nlinsols  = 0;                             # no. of solutions of linear systems
+    ## Print cost statistics if no output argument is given
+    if (nargout == 0)
+      msg = fprintf (1, "Number of successful steps: %d\n", nsteps);
+      msg = fprintf (1, "Number of failed attempts:  %d\n", nfailed);
+      msg = fprintf (1, "Number of function calls:   %d\n", nfevals);
+    endif
+  else
+    havestats = false;
+  endif
+
+  if (nargout == 1)                 # Sort output variables, depends on nargout
+    varargout{1}.x = retvaltime;   # Time stamps are saved in field x
+    varargout{1}.y = retvalresult; # Results are saved in field y
+    varargout{1}.solver = "ode23";  # Solver name is saved in field solver
+    if (haveeventfunction)
+      varargout{1}.ie = event{2};  # Index info which event occured
+      varargout{1}.xe = event{3};  # Time info when an event occured
+      varargout{1}.ye = event{4};  # Results when an event occured
+    endif
+    if (havestats)
+      varargout{1}.stats = struct;
+      varargout{1}.stats.nsteps   = nsteps;
+      varargout{1}.stats.nfailed  = nfailed;
+      varargout{1}.stats.nfevals  = nfevals;
+      varargout{1}.stats.npds     = npds;
+      varargout{1}.stats.ndecomps = ndecomps;
+      varargout{1}.stats.nlinsols = nlinsols;
+    endif
+  elseif (nargout == 2)
+    varargout{1} = retvaltime;     # Time stamps are first output argument
+    varargout{2} = retvalresult;   # Results are second output argument
+  elseif (nargout == 5)
+    varargout{1} = retvaltime;     # Same as (nargout == 2)
+    varargout{2} = retvalresult;   # Same as (nargout == 2)
+    varargout{3} = [];              # LabMat doesn't accept lines like
+    varargout{4} = [];              # varargout{3} = varargout{4} = [];
+    varargout{5} = [];
+    if (haveeventfunction)
+      varargout{3} = event{3};     # Time info when an event occured
+      varargout{4} = event{4};     # Results when an event occured
+      varargout{5} = event{2};     # Index info which event occured
+    endif
+  endif
+endfunction
+
+
+%! ## We are using the "Van der Pol" implementation for all tests that
+%! ## are done for this function. We also define a Jacobian, Events,
+%! ## pseudo-Mass implementation. For further tests we also define a
+%! ## reference solution (computed at high accuracy) and an OutputFcn
+%!function [ydot] = fpol (t, y, varargin) %# The Van der Pol
+%!  ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)];
+%!function [jac] = fjac (t, y, varargin) %# its Jacobian
+%!  jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2];
+%!function [jac] = fjcc (t, y, varargin) %# sparse type
+%!  jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]);
+%!function [val, trm, dir] = feve (t, y, varargin)
+%!  val = fpol (t, y, varargin); # We use the derivatives
+%!  trm = zeros (2,1);             # that's why component 2
+%!  dir = ones (2,1);              # seems to not be exact
+%!function [val, trm, dir] = fevn (t, y, varargin)
+%!  val = fpol (t, y, varargin); # We use the derivatives
+%!  trm = ones (2,1);              # that's why component 2
+%!  dir = ones (2,1);              # seems to not be exact
+%!function [mas] = fmas (t, y)
+%!  mas = [1, 0; 0, 1];            # Dummy mass matrix for tests
+%!function [mas] = fmsa (t, y)
+%!  mas = sparse ([1, 0; 0, 1]);   # A sparse dummy matrix
+%!function [ref] = fref ()         # The computed reference sol
+%!  ref = [0.32331666704577, -1.83297456798624];
+%!function [out] = fout (t, y, flag, varargin)
+%!  if (regexp (char (flag), "init") == 1)
+%!    if (any (size (t) != [2, 1])) error ('\"fout\" step \"init\"'); endif
+%!  elseif (isempty (flag))
+%!    if (any (size (t) != [1, 1])) error ('\"fout\" step \"calc\"'); endif
+%!    out = false;
+%!  elseif (regexp (char (flag), "done") == 1)
+%!    if (any (size (t) != [1, 1])) error ('\"fout\" step \"done\"'); endif
+%!  else error ('\"fout\" invalid flag');
+%!  endif
+%!endfunction
+%!
+%! ## Turn off output of warning messages for all tests, turn them on
+%! ## again if the last test is called
+%!error # input argument number one
+%!  warning ("off", "Octave:InvalidArgument");
+%!  B = ode23 (1, [0 25], [3 15 1]);
+%!error # input argument number two
+%!  B = ode23 (@fpol, 1, [3 15 1]);
+%!error # input argument number three
+%!  B = ode23 (@flor, [0 25], 1);
+%!test # one output argument
+%!  sol = ode23 (@fpol, [0 2], [2 0]);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!  assert (isfield (sol, "solver"));
+%!  assert (sol.solver, "ode23");
+%!test # two output arguments
+%!  [t, y] = ode23 (@fpol, [0 2], [2 0]);
+%!  assert ([t(end), y(end,:)], [2, fref], 1e-3);
+%!test # five output arguments and no Events
+%!  [t, y, xe, ye, ie] = ode23 (@fpol, [0 2], [2 0]);
+%!  assert ([t(end), y(end,:)], [2, fref], 1e-3);
+%!  assert ([ie, xe, ye], []);
+%!test # anonymous function instead of real function
+%!  fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
+%!  sol = ode23 (fvdb, [0 2], [2 0]);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # extra input arguments passed through
+%!  sol = ode23 (@fpol, [0 2], [2 0], 12, 13, "KL");
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # empty Octave structure *but* extra input arguments
+%!  opt = odeset;
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt, 12, 13, "KL");
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!error # strange Octave structure
+%!  opt = struct ("foo", 1);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!test # Solve vdp in fixed step sizes
+%!  sol = ode23 (@fpol, [0:0.1:2], [2 0]);
+%!  assert (sol.x(:), [0:0.1:2]');
+%!  assert (sol.y(end,:), fref, 1e-3);
+%!test # Solve in backward direction starting at t=0
+%!  ref = [-1.205364552835178, 0.951542399860817];
+%!  sol = ode23 (@fpol, [0 -2], [2 0]);
+%!  assert ([sol.x(end), sol.y(end,:)], [-2, ref], 1e-3);
+%!test # Solve in backward direction starting at t=2
+%!  ref = [-1.205364552835178, 0.951542399860817];
+%!  sol = ode23 (@fpol, [2 -2], fref);
+%!  assert ([sol.x(end), sol.y(end,:)], [-2, ref], 1e-3);
+%!test # Solve another anonymous function in backward direction
+%!  ref = [-1, 0.367879437558975];
+%!  sol = ode23 (@(t,y) y, [0 -1], 1);
+%!  assert ([sol.x(end), sol.y(end,:)], ref, 1e-3);
+%!test # Solve another anonymous function below zero
+%!  ref = [0, 14.77810590694212];
+%!  sol = ode23 (@(t,y) y, [-2 0], 2);
+%!  assert ([sol.x(end), sol.y(end,:)], ref, 1e-3);
+%!test # AbsTol option
+%!  opt = odeset ("AbsTol", 1e-5);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # AbsTol and RelTol option
+%!  opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # RelTol and NormControl option -- higher accuracy
+%!  opt = odeset ("RelTol", 1e-8, "NormControl", "on");
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-6);
+%!test # Keeps initial values while integrating
+%!  opt = odeset ("NonNegative", 2);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, 2, 0], 1e-1);
+%!test # Details of OutputSel and Refine can't be tested
+%!  opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!test # Details of OutputSave can't be tested
+%!  opt = odeset ("OutputSave", 1, "OutputSel", 1);
+%!  sla = ode23 (@fpol, [0 2], [2 0], opt);
+%!  opt = odeset ("OutputSave", 2);
+%!  slb = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert (length (sla.x) > length (slb.x))
+%!test # Stats must add further elements in sol
+%!  opt = odeset ("Stats", "on");
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert (isfield (sol, "stats"));
+%!  assert (isfield (sol.stats, "nsteps"));
+%!test # InitialStep option
+%!  opt = odeset ("InitialStep", 1e-8);
+%!  sol = ode23 (@fpol, [0 0.2], [2 0], opt);
+%!  assert ([sol.x(2)-sol.x(1)], [1e-8], 1e-9);
+%!test # MaxStep option
+%!  opt = odeset ("MaxStep", 1e-2);
+%!  sol = ode23 (@fpol, [0 0.2], [2 0], opt);
+%!  assert ([sol.x(5)-sol.x(4)], [1e-2], 1e-2);
+%!test # Events option add further elements in sol
+%!  opt = odeset ("Events", @feve);
+%!  sol = ode23 (@fpol, [0 10], [2 0], opt);
+%!  assert (isfield (sol, "ie"));
+%!  assert (sol.ie, [2; 1; 2; 1]);
+%!  assert (isfield (sol, "xe"));
+%!  assert (isfield (sol, "ye"));
+%!test # Events option, now stop integration
+%!  opt = odeset ("Events", @fevn, "NormControl", "on");
+%!  sol = ode23 (@fpol, [0 10], [2 0], opt);
+%!  assert ([sol.ie, sol.xe, sol.ye], ...
+%!    [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
+%!test # Events option, five output arguments
+%!  opt = odeset ("Events", @fevn, "NormControl", "on");
+%!  [t, y, xe, ye, ie] = ode23 (@fpol, [0 10], [2 0], opt);
+%!  assert ([ie, xe, ye], ...
+%!    [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
+%!test # Jacobian option
+%!  opt = odeset ("Jacobian", @fjac);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # Jacobian option and sparse return value
+%!  opt = odeset ("Jacobian", @fjcc);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.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
+%!  opt = odeset ("Mass", @fmas);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # Mass option as matrix
+%!  opt = odeset ("Mass", eye (2,2));
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # Mass option as sparse matrix
+%!  opt = odeset ("Mass", sparse (eye (2,2)));
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # Mass option as function and sparse matrix
+%!  opt = odeset ("Mass", @fmsa);
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # Mass option as function and MStateDependence
+%!  opt = odeset ("Mass", @fmas, "MStateDependence", "strong");
+%!  sol = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test # Set BDF option to something else than default
+%!  opt = odeset ("BDF", "on");
+%!  [t, y] = ode23 (@fpol, [0 2], [2 0], opt);
+%!  assert ([t(end), y(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", "Octave:InvalidArgument");
+
diff -r 61892a245653 -r b0c8bc80acd4 scripts/ode/ode45.m
--- a/scripts/ode/ode45.m	Sat Dec 12 08:31:22 2015 -0800
+++ b/scripts/ode/ode45.m	Sat Dec 12 18:33:26 2015 +0100
@@ -80,7 +80,7 @@
 ## [@var{t},@var{y}] = ode45 (fvdp, [0 20], [2 0]);
 ## @end group
 ## @end example
-## @seealso{odeset, odeget}
+## @seealso{ode23, odeset, odeget}
 ## @end deftypefn
 
 function varargout = ode45 (fun, trange, init, varargin)
