--- a/scripts/ode/ode45.m	Fri Oct 16 10:31:02 2015 -0700
+++ b/scripts/ode/ode45.m	Sun Oct 18 09:55:41 2015 -0700
@@ -20,388 +20,351 @@
 
 ## -*- 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 (@var{fun}, @var{trange}, @var{init}, @var{ode_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{solution} =} ode45 (@var{fun}, @var{trange}, @var{init}, @dots{})
+## @deftypefnx {Function File} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode45 (@dots{})
+## @deftypefnx {Function File} {@var{solution} =} ode45 (@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{fun} is a function handle, inline function, or string containing the
+## name of the 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
+## evaluated.  Typically, 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}.
+## instances unless the integrate function specified is
+## @command{integrate_n_steps}.
 ##
-## 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}.
+## By default, @code{ode45} uses an adaptive timestep with the
+## @code{integrate_adaptive} algorithm.  The tolerance for the timestep
+## computation may be changed by using the option @qcode{"Tau"}, that has a
+## default value of @math{1e-6}.  If the ODE option @qcode{"TimeStepSize"} is
+## not empty, then the stepper called will be @code{integrate_const}.  If, in
+## addition, the option @qcode{"TimeStepNumber"} is also specified then the
+## integrate function @code{integrate_n_steps} will be used.
 ##
-## If present, the fourth input argument specifies options to the ODE solver.
-## It is a structure typically generated by @code{odeset}.
+## @var{init} contains the initial value for the unknowns.  If it 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}.
 ##
-## The function usually produces just two outputs.  Variable @var{t} is a
+## The optional fourth argument @var{ode_opt} specifies non-default options to
+## the ODE solver.  It is a structure generated by @code{odeset}.
+##
+## The function typically returns 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
+## The output can also be returned as a structure @var{solution} which
+## has field @var{x} containing the time where the solution was evaluated and
+## field @var{y} containing the solution matrix for the times in @var{x}.
+## Use @code{fieldnames (@var{solution})} to see the other fields and additional
+## information returned.
+##
+## If using the @qcode{"Events"} option then three additional outputs may
+## be returned.  @var{te} holds the time when an Event function returned a
+## zero.  @var{ye} holds the value of the solution at time @var{te}.  @var{ie}
+## contains an index indicating which Event function was triggered in the case
+## of multiple Event functions.
+##
+## Example: Solve 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]);
+## fvdp = @@(@var{t},@var{y}) [@var{y}(2); (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)];
+## [@var{t},@var{y}] = ode45 (fvdp, [0 20], [2 0]);
 ## @end group
 ## @end example
 ## @seealso{odeset, odeget}
 ## @end deftypefn
 
-function varargout = ode45 (vfun, vtrange, vinit, varargin)
+function varargout = ode45 (fun, trange, init, varargin)
 
   if (nargin < 3)
     print_usage ();
   endif
-  
-  vorder = 5;  # runge_kutta_45_dorpri uses local extrapolation
-  vsolver = "ode45";
+
+  order = 5;  # runge_kutta_45_dorpri uses local extrapolation
+  solver = "ode45";
 
   if (nargin >= 4)
     if (! isstruct (varargin{1}))
-      ## varargin{1:len} are parameters for vfun
-      vodeoptions = odeset ();
-      vodeoptions.vfunarguments = varargin;
+      ## varargin{1:len} are parameters for fun
+      odeopts = odeset ();
+      odeopts.funarguments = varargin;
     elseif (length (varargin) > 1)
       ## varargin{1} is an ODE options structure vopt
-      vodeoptions = ode_struct_value_check ("ode45", varargin{1}, "ode45");
-      vodeoptions.vfunarguments = {varargin{2:length(varargin)}};
+      odeopts = ode_struct_value_check ("ode45", varargin{1}, "ode45");
+      odeopts.funarguments = {varargin{2:length(varargin)}};
     else  # if (isstruct (varargin{1}))
-      vodeoptions = ode_struct_value_check ("ode45", varargin{1}, "ode45");
-      vodeoptions.vfunarguments = {};
+      odeopts = ode_struct_value_check ("ode45", varargin{1}, "ode45");
+      odeopts.funarguments = {};
     endif
   else  # nargin == 3
-    vodeoptions = odeset (); 
-    vodeoptions.vfunarguments = {};
+    odeopts = odeset ();
+    odeopts.funarguments = {};
   endif
 
-  if (! isvector (vtrange) || ! isnumeric (vtrange))
+  if (! isvector (trange) || ! isnumeric (trange))
     error ("Octave:invalid-input-arg",
-           "second input argument must be a valid vector");
+           "ode45: TRANGE must be a numeric vector");
   endif
 
-  TimeStepNumber = odeget (vodeoptions, "TimeStepNumber", [], "fast");
-  TimeStepSize = odeget (vodeoptions, "TimeStepSize", [], "fast");
-  if (length (vtrange) < 2
+  TimeStepNumber = odeget (odeopts, "TimeStepNumber", [], "fast");
+  TimeStepSize = odeget (odeopts, "TimeStepSize", [], "fast");
+  if (length (trange) < 2
       && (isempty (TimeStepSize) || isempty (TimeStepNumber)))
     error ("Octave:invalid-input-arg",
-           "second input argument must be a valid vector");
-  elseif (vtrange(2) == vtrange(1))
+           "ode45: TRANGE must contain at least 2 elements");
+  elseif (trange(1) == trange(2))
     error ("Octave:invalid-input-arg",
-           "second input argument must be a valid vector");
+           "ode45: invalid time span, TRANGE(1) == TRANGE(2)");
   else
-    vodeoptions.vdirection = sign (vtrange(2) - vtrange(1));
+    odeopts.direction = sign (trange(2) - trange(1));
   endif
-  vtrange = vtrange(:);
+  trange = trange(:);
 
-  if (! isvector (vinit) || ! isnumeric (vinit))
+  if (! isvector (init) || ! isnumeric (init))
     error ("Octave:invalid-input-arg",
-           "third input argument must be a valid numerical value");
+           "ode45: INIT must be a numeric vector");
   endif
-  vinit = vinit(:);
+  init = init(:);
 
-  if (ischar (vfun))
+  if (ischar (fun))
     try
-      vfun = str2func (vfun);
+      fun = str2func (fun);
     catch
       warning (lasterr);
     end_try_catch
   endif
-  if (! isa (vfun, "function_handle"))
+  if (! isa (fun, "function_handle"))
     error ("Octave:invalid-input-arg",
-           "first input argument must be a valid function handle");
+           "ode45: FUN must be a valid function handle");
   endif
 
-  ## Start preprocessing, have a look which options are set in vodeoptions,
+  ## Start preprocessing, have a look which options are set in odeopts,
   ## check if an invalid or unused option is set
   if (isempty (TimeStepNumber) && isempty (TimeStepSize))
     integrate_func = "adaptive";
-    vodeoptions.vstepsizefixed = false;
+    odeopts.stepsizefixed = false;
   elseif (! isempty (TimeStepNumber) && ! isempty (TimeStepSize))
     integrate_func = "n_steps";
-    vodeoptions.vstepsizefixed = true;
-    if (sign (TimeStepSize) != vodeoptions.vdirection)
+    odeopts.stepsizefixed = true;
+    if (sign (TimeStepSize) != odeopts.direction)
       warning ("Octave:invalid-input-arg",
-               "option 'TimeStepSize' has a wrong sign",
-               "it will be corrected automatically");
+               ["ode45: option 'TimeStepSize' has the wrong sign, ", ...
+                "but will be corrected automatically\n"]);
       TimeStepSize = -TimeStepSize;
     endif
   elseif (isempty (TimeStepNumber) && ! isempty (TimeStepSize))
     integrate_func = "const";
-    vodeoptions.vstepsizefixed = true;
-    if (sign (TimeStepSize) != vodeoptions.vdirection)
+    odeopts.stepsizefixed = true;
+    if (sign (TimeStepSize) != odeopts.direction)
       warning ("Octave:invalid-input-arg",
-               "option 'TimeStepSize' has a wrong sign",
-               "it will be corrected automatically");
+               ["ode45: option 'TimeStepSize' has the wrong sign, ",
+                "but will be corrected automatically\n"]);
       TimeStepSize = -TimeStepSize;
     endif
   else
     warning ("Octave:invalid-input-arg",
-             "assuming an adaptive integrate function");
+             "ode45: assuming an adaptive integrate function\n");
     integrate_func = "adaptive";
   endif
-  vodeoptions.TimeStepSize = TimeStepSize;
-  vodeoptions.TimeStepNumber = TimeStepNumber;
-
-  ## 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 ("Octave:invalid-input-arg",
-             "Option 'RelTol' not set, new value %f is used",
-             vodeoptions.RelTol);
-  elseif (! isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed)
-    warning ("Octave:invalid-input-arg",
-             "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;
+  if (isempty (odeopts.RelTol) && ! odeopts.stepsizefixed)
+    odeopts.RelTol = 1e-3;
     warning ("Octave:invalid-input-arg",
-             "Option 'AbsTol' not set, new value %f is used",
-             vodeoptions.AbsTol);
-  elseif (! isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed)
+             "ode45: option 'RelTol' not set, new value %f will be used\n",
+             odeopts.RelTol);
+  elseif (! isempty (odeopts.RelTol) && odeopts.stepsizefixed)
     warning ("Octave:invalid-input-arg",
-             "Option 'AbsTol' will be ignored if fixed time stamps are given");
-  else
-    vodeoptions.AbsTol = vodeoptions.AbsTol(:); # Create column vector
+             ["ode45: option 'RelTol' is ignored", ...
+              " when fixed time stamps are given\n"]);
   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 ("Octave:invalid-input-arg",
-               "Option 'NonNegative' will be ignored if mass matrix is set");
-    endif
-  else 
-    vodeoptions.vhavenonnegative = false;
+  if (isempty (odeopts.AbsTol) && ! odeopts.stepsizefixed)
+    odeopts.AbsTol = 1e-6;
+    warning ("Octave:invalid-input-arg",
+             "ode45: option 'AbsTol' not set, new value %f will be used\n",
+             odeopts.AbsTol);
+  elseif (! isempty (odeopts.AbsTol) && odeopts.stepsizefixed)
+    warning ("Octave:invalid-input-arg",
+             ["ode45: option 'AbsTol' is ignored", ...
+              " when fixed time stamps are given\n"]);
+  else
+    odeopts.AbsTol = odeopts.AbsTol(:);  # Create column vector
   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;
+  odeopts.normcontrol = strcmp (odeopts.NormControl, "on");
+
+  if (! isempty (odeopts.NonNegative))
+    if (isempty (odeopts.Mass))
+      odeopts.havenonnegative = true;
+    else
+      odeopts.havenonnegative = false;
+      warning ("Octave:invalid-input-arg",
+               ["ode45: option 'NonNegative' is ignored", ...
+                " when mass matrix is set\n"]);
+    endif
+  else
+    odeopts.havenonnegative = false;
   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 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;
+  if (isempty (odeopts.OutputFcn) && nargout == 0)
+    odeopts.OutputFcn = @odeplot;
+    odeopts.haveoutputfunction = true;
+  else
+    odeopts.haveoutputfunction = ! isempty (odeopts.OutputFcn);
   endif
 
-  ## Implementation of the option Stats has been finished.
-  ## This option can be set by the user to another value than default value.
+  odeopts.haveoutputselection = ! isempty (odeopts.OutputSel);
 
-  ## 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 ("Octave:invalid-input-arg",
-             "option 'InitialStep' not set, estimated value %f is used",
-             vodeoptions.InitialStep);
-  elseif (isempty (vodeoptions.InitialStep))
-    vodeoptions.InitialStep = TimeStepSize;
+  if (odeopts.Refine > 0)
+    odeopts.haverefine = true;
+  else
+    odeopts.haverefine = false;
   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;
+  if (isempty (odeopts.InitialStep) && strcmp (integrate_func, "adaptive"))
+    odeopts.InitialStep = ...
+      odeopts.direction * starting_stepsize (order, fun, trange(1),
+                                                  init,
+                                                  odeopts.AbsTol,
+                                                  odeopts.RelTol,
+                                                  odeopts.normcontrol);
     warning ("Octave:invalid-input-arg",
-             "Option 'MaxStep' not set, new value %f is used",
-             vodeoptions.MaxStep);
+             ["ode45: option 'InitialStep' not set,", ...
+              " estimated value %f will be used\n"],
+             odeopts.InitialStep);
+  elseif (isempty (odeopts.InitialStep))
+    odeopts.InitialStep = TimeStepSize;
   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;
+  if (isempty (odeopts.MaxStep) && ! odeopts.stepsizefixed)
+    odeopts.MaxStep = abs (trange(end) - trange(1)) / 10;
+    warning ("Octave:invalid-input-arg",
+             "ode45: option 'MaxStep' not set, new value %f will be used\n",
+             odeopts.MaxStep);
   endif
 
+  odeopts.haveeventfunction = ! isempty (odeopts.Events);
+
   ## 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))
+  if (! isempty (odeopts.Jacobian))
     warning ("Octave:invalid-input-arg",
-             "option 'Jacobian' will be ignored by this solver");
+             "ode45: option 'Jacobian' is ignored by this solver\n");
   endif
 
-  if (! isequal (vodeoptions.JPattern, vodetemp.JPattern))
+  if (! isempty (odeopts.JPattern))
     warning ("Octave:invalid-input-arg",
-             "option 'JPattern' will be ignored by this solver");
+             "ode45: option 'JPattern' is ignored by this solver\n");
   endif
 
-  if (! isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
+  if (! isempty (odeopts.Vectorized))
     warning ("Octave:invalid-input-arg",
-             "option 'Vectorized' will be ignored by this solver");
+             "ode45: option 'Vectorized' is ignored by this solver\n");
   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
+  if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass))
+    havemasshandle = false;
+    mass = odeopts.Mass;  # constant mass
+  elseif (isa (odeopts.Mass, "function_handle"))
+    havemasshandle = 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);
+    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 (vodeoptions.MStateDependence, "none"))
-    vmassdependence = false;
-  else
-    vmassdependence = true;
+  massdependence = ! strcmp (odeopts.MStateDependence, "none");
+
+  ## Other options that are not used by this solver.
+  if (! isempty (odeopts.MvPattern))
+    warning ("Octave:invalid-input-arg",
+             "ode45: option 'MvPattern' is ignored by this solver\n");
   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))
+  if (! isempty (odeopts.MassSingular))
     warning ("Octave:invalid-input-arg",
-             "option 'MvPattern' will be ignored by this solver");
+             "ode45: option 'MassSingular' is ignored by this solver\n");
   endif
 
-  if (! isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
+  if (! isempty (odeopts.InitialSlope))
     warning ("Octave:invalid-input-arg",
-             "option 'MassSingular' will be ignored by this solver");
+             "ode45: option 'InitialSlope' is ignored by this solver\n");
   endif
 
-  if (! isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
+  if (! isempty (odeopts.MaxOrder))
     warning ("Octave:invalid-input-arg",
-             "option 'InitialSlope' will be ignored by this solver");
+             "ode45: option 'MaxOrder' is ignored by this solver\n");
   endif
 
-  if (! isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
+  if (! isempty (odeopts.BDF))
     warning ("Octave:invalid-input-arg",
-             "option 'MaxOrder' will be ignored by this solver");
-  endif
-
-  if (! isequal (vodeoptions.BDF, vodetemp.BDF))
-    warning ("Octave:invalid-input-arg",
-             "option 'BDF' will be ignored by this solver");
+             "ode45: option 'BDF' is ignored by this solver\n");
   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{:});
+
+  if (havemasshandle)   # Handle only the dynamic mass matrix,
+    if (massdependence) # constant mass matrices have already
+      mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:});
+      fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ...
+             \ fun (t, x, odeopts.funarguments{:});
+    else                 # if (massdependence == false)
+      mass = @(t) odeopts.Mass (t, odeopts.funarguments{:});
+      fun = @(t,x) mass (t, odeopts.funarguments{:}) ...
+             \ fun (t, x, odeopts.funarguments{:});
     endif
   endif
 
   switch (integrate_func)
     case "adaptive"
       solution = integrate_adaptive (@runge_kutta_45_dorpri,
-                                     vorder, vfun, vtrange, vinit, SubOpts);
+                                     order, fun, trange, init, odeopts);
     case "n_steps"
       solution = integrate_n_steps (@runge_kutta_45_dorpri,
-                                    vfun, vtrange(1), vinit,
-                                    TimeStepSize, TimeStepNumber, SubOpts);
+                                    fun, trange(1), init,
+                                    TimeStepSize, TimeStepNumber, odeopts);
     case "const"
       solution = integrate_const (@runge_kutta_45_dorpri,
-                                  vfun, vtrange, vinit,
-                                  TimeStepSize, SubOpts);
+                                  fun, trange, init,
+                                  TimeStepSize, odeopts);
   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{:});
+  if (odeopts.haveoutputfunction)  # Cleanup plotter
+    feval (odeopts.OutputFcn, solution.t(end),
+           solution.x(end,:)', "done", odeopts.funarguments{:});
   endif
-  if (vodeoptions.vhaveeventfunction)   # Cleanup event function handling
-    odepkg_event_handle (vodeoptions.Events, solution.t(end),
+  if (odeopts.haveeventfunction)   # Cleanup event function handling
+    ode_event_handler (odeopts.Events, solution.t(end),
                          solution.x(end,:)', "done",
-                         vodeoptions.vfunarguments{:});
+                         odeopts.funarguments{:});
   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   = 6 * (solution.vcntcycles-1) + 1;     # number of ode evaluations
-    vndecomps  = 0;                                   # number of LU decompositions
-    vnpds      = 0;                                   # number of partial derivatives
-    vnlinsols  = 0;                                   # no. of linear systems solutions
+  if (strcmp (odeopts.Stats, "on"))
+    havestats = true;
+    nsteps    = solution.cntloop-2;                 # cntloop from 2..end
+    nfailed   = (solution.cntcycles-1)-(nsteps)+1;  # cntcycl from 1..end
+    nfevals   = 6 * (solution.cntcycles-1) + 1;     # number of ode evaluations
+    ndecomps  = 0;  # number of LU decompositions
+    npds      = 0;  # number of partial derivatives
+    nlinsols  = 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);
+      printf ("Number of successful steps: %d\n", nsteps);
+      printf ("Number of failed attempts:  %d\n", nfailed);
+      printf ("Number of function calls:   %d\n", nfevals);
     endif
   else
-    vhavestats = false;
+    havestats = false;
   endif
 
   if (nargout == 2)
@@ -410,29 +373,29 @@
   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
+    varargout{1}.solver = solver;   # Solver name is saved in field solver
+    if (odeopts.haveeventfunction)
+      varargout{1}.ie = solution.event{2};  # Index info which event occurred
+      varargout{1}.xe = solution.event{3};  # Time info when an event occurred
+      varargout{1}.ye = solution.event{4};  # Results when an event occurred
     endif
-    if (vhavestats)
+    if (havestats)
       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;
+      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 == 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
+    if (odeopts.haveeventfunction)
+      varargout{3} = solution.event{3};  # Time info when an event occurred
+      varargout{4} = solution.event{4};  # Results when an event occurred
+      varargout{5} = solution.event{2};  # Index info which event occurred
     endif
   endif
 
@@ -443,210 +406,211 @@
 ## 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)];
+%!function ydot = fpol (t, y)  # The Van der Pol
+%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)];
 %!endfunction
-%!function [vref] = fref ()        # The computed reference solution
-%! vref = [0.32331666704577, -1.83297456798624];
+%!function ref = fref ()         # The computed reference solution
+%! ref = [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];
+%!function jac = fjac (t, y, varargin) # its Jacobian
+%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(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])
+%!function jac = fjcc (t, y, varargin) # sparse type
+%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(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
+%!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
 %!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
+%!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
 %!endfunction
-%!function [vmas] = fmas (vt, vy, varargin)
-%! vmas = [1, 0; 0, 1];            # Dummy mass matrix for tests
+%!function mas = fmas (t, y, varargin)
+%! mas = [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
+%!function mas = fmsa (t, y, varargin)
+%! mas = 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
+%!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 vflag');
+%!   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  # ouput argument
-%! warning ("off", "Octave:invalid-input-arg");
-%! B = ode45 (1, [0 25], [3 15 1]);
-%!error  # input argument number one
-%! [vt, vy] = ode45 (1, [0 25], [3 15 1]);
+%! ## Turn off output of warning messages for all tests,
+%! ## turn them on again after the last test is called
+%!test
+%! warning ("off", "Octave:invalid-input-arg", "local");
+%!error  # first input arg is not a function
+%! [t, y] = ode45 (1, [0 25], [3 15 1]);
 %!error  # input argument number one as name of non existing function
-%! [vt, vy] = ode45 ("non-existing-function", [0 25], [3 15 1]);
+%! [t, y] = ode45 ("non-existing-function", [0 25], [3 15 1]);
 %!error  # input argument number two
-%! [vt, vy] = ode45 (@fpol, 1, [3 15 1]);
+%! [t, y] = 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);
+%! [t, y] = ode45 (@fpol, [0 2], [2 0]);
+%! assert ([t(end), y(end,:)], [2, fref], 1e-2);
 %!test  # not too many steps
-%! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
-%! assert (size (vt) < 20);
+%! [t, y] = ode45 (@fpol, [0 2], [2 0]);
+%! assert (size (t) < 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);
+%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
+%! [t, y] = ode45 (fvdb, [0 2], [2 0]);
+%! assert ([t(end), y(end,:)], [2, fref], 1e-2);
 %!test  # string instead of function
-%! [vt, vy] = ode45 ("fpol", [0 2], [2 0]);
-%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
+%! [t, y] = ode45 ("fpol", [0 2], [2 0]);
+%! assert ([t(end), y(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);
+%! [t, y] = ode45 (@fpol, [0 2], [2 0], 12, 13, "KL");
+%! assert ([t(end), y(end,:)], [2, fref], 1e-2);
 %!test  # empty ODEOPT 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);
+%! opt = odeset;
+%! [t, y] = ode45 (@fpol, [0 2], [2 0], opt, 12, 13, "KL");
+%! assert ([t(end), y(end,:)], [2, fref], 1e-2);
 %!error  # strange ODEOPT structure
-%! vopt = struct ("foo", 1);
-%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
+%! opt = struct ("foo", 1);
+%! [t, y] = ode45 (@fpol, [0 2], [2 0], opt);
 %!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);
+%! opt = odeset("TimeStepSize", 0.1);
+%! [t, y] = ode45 (@fpol, [0,2], [2 0], opt);
+%! assert (t(:), [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);
+%! [t, y] = ode45 (@(t,y) y, [-2 0], 2);
+%! assert ([t(end), y(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);
+%! opt = odeset ("InitialStep", 1e-8);
+%! [t, y] = ode45 (@fpol, [0 0.2], [2 0], opt);
+%! assert ([t(2)-t(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);
+%! opt = odeset ("MaxStep", 1e-3);
+%! sol = ode45 (@fpol, [0 0.2], [2 0], opt);
+%! assert ([sol.x(5)-sol.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);
+%! sol = ode45 (@fpol, [0 1 2], [2 0]);
+%! assert (any((sol.x-1) == 0));
+%! assert ([sol.x(end), sol.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);
+%! sol = ode45 (@fpol, [0 -2], [2 0]);
+%! assert ([sol.x(end), sol.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);
+%! sol = ode45 (@fpol, [2 -2], fref);
+%! assert ([sol.x(end), sol.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);
+%! sol = ode45 (@fpol, [2 0 -2], fref);
+%! idx = find(sol.x < 0, 1, "first") - 1;
+%! assert ([sol.x(idx), sol.y(idx,:)], [0 2 0], 1e-2);
+%! assert ([sol.x(end), sol.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);
+%! sol = ode45 (@(t,y) y, [0 -1], 1);
+%! assert ([sol.x(end), sol.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);
+%! sol = ode45 (@(t,y) y, [-2 0], 2);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("MaxStep", 1e-3);
+%! sol = ode45 (@fpol, [0 -2], [2 0], opt);
+%! assert ([abs(sol.x(8)-sol.x(7))], [1e-3], 1e-3);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("AbsTol", 1e-5);
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8);
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("RelTol", 1e-8, "NormControl", "on");
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("NonNegative", 2);
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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  # 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"));
+%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%!test  # Stats must add further elements in sol
+%! opt = odeset ("Stats", "on");
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert (isfield (sol, "stats"));
+%! assert (isfield (sol.stats, "nsteps"));
+%!test  # Events option add further elements in sol
+%! opt = odeset ("Events", @feve);
+%! sol = ode45 (@fpol, [0 10], [2 0], opt);
+%! assert (isfield (sol, "ie"));
+%! assert (sol.ie(1), 2);
+%! assert (isfield (sol, "xe"));
+%! assert (isfield (sol, "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], 
+%! opt = odeset ("Events", @fevn, "NormControl", "on");
+%! sol = ode45 (@fpol, [0 10], [2 0], opt);
+%! assert ([sol.ie, sol.xe, sol.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],
+%! opt = odeset ("Events", @fevn, "NormControl", "on");
+%! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt);
+%! assert ([vie, vxe, ye],
 %!         [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);
+%! opt = odeset ("Jacobian", @fjac);
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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
-%!
+%! opt = odeset ("Jacobian", @fjcc);
+%! sol = ode45 (@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  # 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);
+%! opt = odeset ("Mass", @fmas);
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("Mass", eye (2,2));
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("Mass", sparse (eye (2,2)));
+%! sol = ode45 (@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
-%! vopt = odeset ("Mass", @fmsa);
-%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
+%! opt = odeset ("Mass", @fmsa);
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.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);
+%! opt = odeset ("Mass", @fmas, "MStateDependence", "strong");
+%! sol = ode45 (@fpol, [0 2], [2 0], opt);
+%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
+%!test  # Set BDF option to something other than default
+%! opt = odeset ("BDF", "on");
+%! [t, y] = ode45 (@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
+
 %!test
-%!## test for MvPattern option is missing
-%!## test for InitialSlope option is missing
-%!## test for MaxOrder option is missing
 %!
-%! warning ("on", "Octave:invalid-input-arg");
+%! #warning ("on", "Octave:invalid-input-arg");