changeset 20581:e368ce72a844

maint: Use Octave coding conventions for ode* functions. Use '!' instead of '~' for logical not. Use '##' only for comments which standalone on a line. Use double quotes instead of single quotes where possible. Wrap long lines to < 80 characters. * ode45.m, odeget.m, odeset.m: Use Octave coding conventions.
author Rik <rik@octave.org>
date Sat, 03 Oct 2015 22:10:45 -0700
parents 25623ef2ff4f
children 235d059cf329
files scripts/ode/ode45.m scripts/ode/odeget.m scripts/ode/odeset.m
diffstat 3 files changed, 403 insertions(+), 397 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ode/ode45.m	Sat Oct 03 21:03:16 2015 -0700
+++ b/scripts/ode/ode45.m	Sat Oct 03 22:10:45 2015 -0700
@@ -72,12 +72,12 @@
 ## @seealso{odeset, odeget}
 ## @end deftypefn
 
-function [varargout] = ode45 (vfun, vslot, vinit, varargin)
+function varargout = ode45 (vfun, vslot, vinit, varargin)
 
-  vorder = 5; % runge_kutta_45_dorpri uses local extrapolation
+  vorder = 5; # runge_kutta_45_dorpri uses local extrapolation
   vsolver = "ode45";
 
-  if (nargin == 0) ## Check number and types of all input arguments
+  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");
@@ -88,7 +88,7 @@
   endif
   
   if (nargin >= 4)
-    if (~isstruct (varargin{1}))
+    if (! isstruct (varargin{1}))
       ## varargin{1:len} are parameters for vfun
       vodeoptions = odeset;
       vodeoptions.vfunarguments = varargin;
@@ -96,16 +96,16 @@
       ## 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}))
+    else  # if (isstruct (varargin{1}))
       vodeoptions = odepkg_structure_check (varargin{1}, "ode45");
       vodeoptions.vfunarguments = {};
     endif
-  else ## if (nargin == 3)
+  else  # nargin == 3
     vodeoptions = odeset; 
     vodeoptions.vfunarguments = {};
   endif
 
-  if (~isvector (vslot) || ~isnumeric (vslot))
+  if (! isvector (vslot) || ! isnumeric (vslot))
     error ("OdePkg:InvalidArgument", ...
            "second input argument must be a valid vector");
   endif
@@ -123,39 +123,40 @@
   endif
   vslot = vslot(:);
 
-  if (~isvector (vinit) || ~isnumeric (vinit))
+  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") || isa (vfun, "inline"))
+  if (! (isa (vfun, "function_handle") || isa (vfun, "inline")))
     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
+  ## 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))
+  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", ...
+               "option 'TimeStepSize' has a wrong sign", ...
                "it will be corrected automatically");
       vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
     endif
-  elseif (isempty (vodeoptions.TimeStepNumber) && ~isempty (vodeoptions.TimeStepSize))
+  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", ...
+               "option 'TimeStepSize' has a wrong sign", ...
                "it will be corrected automatically");
       vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
     endif
@@ -168,43 +169,43 @@
   ## 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)
+  ## 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", ...
+             "Option 'RelTol' not set, new value %f is used", ...
              vodeoptions.RelTol);
-  elseif (~isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed)
+  elseif (! isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed)
     warning ("OdePkg:InvalidArgument", ...
-             "Option ''RelTol'' will be ignored if fixed time stamps are given");
+             "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)
+  ## 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", ...
+             "Option 'AbsTol' not set, new value %f is used", ...
              vodeoptions.AbsTol);
-  elseif (~isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed)
+  elseif (! isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed)
     warning ("OdePkg:InvalidArgument", ...
-             "Option ''AbsTol'' will be ignored if fixed time stamps are given");
+             "Option 'AbsTol' will be ignored if fixed time stamps are given");
   else
-    vodeoptions.AbsTol = vodeoptions.AbsTol(:); ## Create column vector
+    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.
+  ## 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))
+  ## 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
@@ -216,8 +217,8 @@
     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.
+  ## 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;
@@ -227,39 +228,41 @@
     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))
+  ## 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.
+  ## 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.
+  ## 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 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.
+  ## 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, vslot(1), vinit, vodeoptions.AbsTol, ...
-                                                 vodeoptions.RelTol, vodeoptions.vnormcontrol);
+    vodeoptions.InitialStep = ...
+      vodeoptions.vdirection * starting_stepsize (vorder, vfun, vslot(1), vinit,
+                                                  vodeoptions.AbsTol,
+                                                  vodeoptions.RelTol,
+                                                  vodeoptions.vnormcontrol);
     warning ("OdePkg:InvalidArgument", ...
-             "option ''InitialStep'' not set, estimated value %f is used", ...
+             "option 'InitialStep' not set, estimated value %f is used", ...
              vodeoptions.InitialStep);
   elseif(isempty (vodeoptions.InitialStep))
     vodeoptions.InitialStep = odeget (vodeoptions, "TimeStepSize");
@@ -267,105 +270,104 @@
 
   ## 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)
+  if (isempty (vodeoptions.MaxStep) && ! vodeoptions.vstepsizefixed)
     vodeoptions.MaxStep = abs (vslot(end) - vslot(1)) / 10;
     warning ("OdePkg:InvalidArgument", ...
-             "Option ''MaxStep'' not set, new value %f is used", ...
+             "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))
+  ## 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))
+  ## 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");
+             "option 'Jacobian' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
+  if (! isequal (vodeoptions.JPattern, vodetemp.JPattern))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''JPattern'' will be ignored by this solver");
+             "option 'JPattern' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
+  if (! isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''Vectorized'' will be ignored by this solver");
+             "option 'Vectorized' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
+  if (! isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''NewtonTol'' will be ignored by this solver");
+             "option 'NewtonTol' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.MaxNewtonIterations,vodetemp.MaxNewtonIterations))
+  if (! isequal (vodeoptions.MaxNewtonIterations,vodetemp.MaxNewtonIterations))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''MaxNewtonIterations'' will be ignored by this solver");
+             "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))
+  ## 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
+    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);
+    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.
+  ## 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(s) is/are ignored.
-  if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
+  ## 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");
+             "option 'MvPattern' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
+  if (! isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''MassSingular'' will be ignored by this solver");
+             "option 'MassSingular' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
+  if (! isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''InitialSlope'' will be ignored by this solver");
+             "option 'InitialSlope' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
+  if (! isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''MaxOrder'' will be ignored by this solver");
+             "option 'MaxOrder' will be ignored by this solver");
   endif
 
-  if (~isequal (vodeoptions.BDF, vodetemp.BDF))
+  if (! isequal (vodeoptions.BDF, vodetemp.BDF))
     warning ("OdePkg:InvalidArgument", ...
-             "option ''BDF'' will be ignored by this solver");
+             "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
+  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)
+    else                 # if (vmassdependence == false)
       vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:});
       vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ...
              \ vfun (t, x, vodeoptions.vfunarguments{:});
@@ -388,24 +390,25 @@
   endswitch
 
   ## Postprocessing, do whatever when terminating integration algorithm
-  if (vodeoptions.vhaveoutputfunction) ## Cleanup plotter
+  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
+  if (vodeoptions.vhaveeventfunction)   # Cleanup event function handling
     odepkg_event_handle (vodeoptions.Events, solution.t(end), ...
-                         solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
+                         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)-(solution.vcntloop-2)+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 solutions of linear systems
+    vnsteps    = solution.vcntloop-2;        # vcntloop from 2..end
+    vnfailed   = (solution.vcntcycles-1)-(solution.vcntloop-2)+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)
       vmsg = fprintf (1, "Number of successful steps: %d\n", vnsteps);
@@ -416,14 +419,14 @@
     vhavestats = false;
   endif
 
-  if (nargout == 1)                 ## Sort output variables, depends on nargout
-    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 (nargout == 1)                 # Sort output variables, depends on nargout
+    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 occured
-      varargout{1}.xe = solution.vevent{3};  ## Time info when an event occured
-      varargout{1}.ye = solution.vevent{4};  ## Results when an event occured
+      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;
@@ -435,228 +438,235 @@
       varargout{1}.stats.nlinsols = vnlinsols;
     endif
   elseif (nargout == 2)
-    varargout{1} = solution.t;     ## Time stamps are first output argument
-    varargout{2} = solution.x;   ## Results are second output argument
+    varargout{1} = solution.t;   # Time stamps are first output argument
+    varargout{2} = solution.x;   # Results are second output argument
   elseif (nargout == 5)
-    varargout{1} = solution.t;     ## Same as (nargout == 2)
-    varargout{2} = solution.x;   ## Same as (nargout == 2)
-    varargout{3} = [];              ## LabMat doesn't accept lines like
-    varargout{4} = [];              ## varargout{3} = varargout{4} = [];
+    varargout{1} = solution.t;   # Same as (nargout == 2)
+    varargout{2} = solution.x;   # Same as (nargout == 2)
+    varargout{3} = [];           # LabMat doesn't accept lines like
+    varargout{4} = [];           # varargout{3} = varargout{4} = [];
     varargout{5} = [];
     if (vodeoptions.vhaveeventfunction) 
-      varargout{3} = solution.vevent{3};     ## Time info when an event occured
-      varargout{4} = solution.vevent{4};     ## Results when an event occured
-      varargout{5} = solution.vevent{2};     ## Index info which event occured
+      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)];
+
+## 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 sol
-%!  vref = [0.32331666704577, -1.83297456798624];
+%!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];
-%!function [vjac] = fjcc (vt, vy, varargin) ## sparse type
-%!  vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2])
+%!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
+%! 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
+%! 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
+%! 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
+%! 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"'); end
-%!  elseif (isempty (vflag))
-%!    if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
-%!    vout = false;
-%!  elseif (regexp (char (vflag), 'done') == 1)
-%!    if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
-%!  else error ('"fout" invalid vflag');
-%!  end
+%! 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);
+%!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 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  # 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
+%!## test for MvPattern option is missing
+%!## test for InitialSlope option is missing
+%!## test for MaxOrder option is missing
 %!
-%!  warning ('on', 'OdePkg:InvalidArgument');
+%! warning ("on", "OdePkg:InvalidArgument");
 
-## Local Variables: ***
-## mode: octave ***
-## End: ***
--- a/scripts/ode/odeget.m	Sat Oct 03 21:03:16 2015 -0700
+++ b/scripts/ode/odeget.m	Sat Oct 03 22:10:45 2015 -0700
@@ -28,9 +28,11 @@
 ## argument @var{ode_opt} is of type structure array and the second input
 ## argument @var{field} is of type string then return the option value
 ## @var{res} that is specified by the option name @var{field} in the ODE
-## option structure @var{ode_opt}.  Optionally if this function is called
-## with a third input argument then return the default value @var{default} if
-## @var{field} is not set in the structure @var{ode_opt}.
+## option structure @var{ode_opt}.
+##
+## Optionally if this function is called with a third input argument then
+## return the default value @var{default} if @var{field} is not set in the
+## structure @var{ode_opt}.
 ## @seealso{odeset}
 ## @end deftypefn
 
@@ -111,10 +113,10 @@
     if (size (pos, 1) == 0) # no matching for the given option
       if (nargin == 2)
         error ("OdePkg:InvalidArgument",
-               "invalid property. No property found with name ''%s''.", field);
+               "invalid property. No property found with name '%s'.", field);
       endif
       warning ("odeget:NoExactMatching",
-               "no property found with name ''%s''. ",
+               "no property found with name '%s'. ",
                "Assuming default value.", field);
       res = default;
       return
@@ -124,8 +126,8 @@
       if (! strcmp (lower (deblank (field)),
                     lower (deblank (options(pos,:)))) )
         warning ("OdePkg:InvalidArgument",
-                 "no exact matching for ''%s''. ",
-                 "Assuming you was intending ''%s''.",
+                 "no exact matching for '%s'. ",
+                 "Assuming you was intending '%s'.",
                  field, deblank (options(pos,:)));
       endif
       res = ode_opt.(deblank (options(pos,:)));
@@ -138,7 +140,7 @@
 
     ## if there are more matching, ask the user to be more precise
     warning ("OdePkg:InvalidArgument", ...
-             "no exact matching for ''%s''. %d possible fields were found.",
+             "no exact matching for '%s'. %d possible fields were found.",
              field, size (pos, 1));
     for j = 1:(size (pos, 1))
       fprintf ("%s\n", deblank (options(pos(j),:)));
@@ -151,27 +153,25 @@
 
 endfunction
 
-%! ## Turn off output of warning messages for all tests, turn them on
-%! ## again if the last test is called
-%!  warning ('off', 'OdePkg:InvalidArgument');
-%!test assert (odeget (odeset (), 'RelTol'), []);
-%!test assert (odeget (odeset (), 'RelTol', 10), 10);
-%!test assert (odeget (odeset (), 'Stats'), []);
-%!test assert (odeget (odeset (), 'Stats', 'on'), 'on');
-%!test assert (odeget (odeset (), 'AbsTol', 1.e-6, 'fast'), []);
-%!test assert (odeget (odeset (), 'AbsTol', 1.e-6, 'fast_not_empty'), 1.e-6);
-%!test assert (odeget (odeset (), 'AbsTol', 1e-9), 1e-9);
-%!
-%!  warning ('on', 'OdePkg:InvalidArgument');
 
 %!demo
 %! # Return the manually changed value RelTol of the OdePkg options
-%! # strutcure A. If RelTol wouldn't have been changed then an
+%! # strutcure A.  If RelTol wouldn't have been changed then an
 %! # empty matrix value would have been returned.
 %!
-%! A = odeset ('RelTol', 1e-1, 'AbsTol', 1e-2);
-%! odeget (A, 'RelTol', [])
+%! A = odeset ("RelTol", 1e-1, "AbsTol", 1e-2);
+%! odeget (A, "RelTol", [])
 
-## Local Variables: ***
-## mode: octave ***
-## End: ***
+%! ## Turn off output of warning messages for all tests, turn them on
+%! ## again if the last test is called
+%!  warning ("off", "OdePkg:InvalidArgument");
+%!test assert (odeget (odeset (), "RelTol"), []);
+%!test assert (odeget (odeset (), "RelTol", 10), 10);
+%!test assert (odeget (odeset (), "Stats"), []);
+%!test assert (odeget (odeset (), "Stats", "on"), "on");
+%!test assert (odeget (odeset (), "AbsTol", 1.e-6, "fast"), []);
+%!test assert (odeget (odeset (), "AbsTol", 1.e-6, "fast_not_empty"), 1.e-6);
+%!test assert (odeget (odeset (), "AbsTol", 1e-9), 1e-9);
+%!
+%!  warning ("on", "OdePkg:InvalidArgument");
+
--- a/scripts/ode/odeset.m	Sat Oct 03 21:03:16 2015 -0700
+++ b/scripts/ode/odeset.m	Sat Oct 03 22:10:45 2015 -0700
@@ -17,7 +17,6 @@
 ## along with Octave; see the file COPYING.  If not, see
 ## <http://www.gnu.org/licenses/>.
 
-
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} odeset ()
 ## @deftypefnx {Function File} {@var{odestruct} =} odeset (@var{"field1"}, @var{value1}, @var{"field2"}, @var{value2}, @dots{})
@@ -53,10 +52,9 @@
 function opt = odeset (varargin)
 
   ## Check number and types of all input arguments
-  if ((nargin == 0)
-      && (nargout == 0))
+  if (nargin == 0 && nargout == 0)
     print_options;
-    return
+    return;
   endif
 
   ## Creating a vector of OdePkg possible fields
@@ -81,9 +79,8 @@
   opt.Refine = 0;
   opt.OutputSave = 1;
 
-  if ((nargin == 0)
-      && (nargout == 1))
-    return
+  if (nargin == 0 && nargout == 1)
+    return;
   endif
 
   ode_fields = fieldnames (opt);
@@ -102,14 +99,14 @@
         pos = fuzzy_compare (name, fields);
         if (size (pos, 1) == 0)
           warning ("OdePkg:InvalidArgument",
-                   "no property found with name ''%s''", name);
+                   "no property found with name '%s'", name);
         endif
 
         if (size (pos, 1) == 1)
           if (! strcmp (lower (deblank (name)),
                         lower (deblank (fields(pos,:)))))
             warning ("OdePkg:InvalidArgument", "no exact matching for ",
-                     "''%s''. Assuming you were intending ''%s''",
+                     "'%s'. Assuming you were intending '%s'",
                      name, deblank (fields(pos,:)));
           endif
 
@@ -119,7 +116,7 @@
 
         ## if there are more matching, ask the user to be more precise
         warning ("OdePkg:InvalidArgument",
-                 "no exact matching for ''%s''. %d possible fields were found",
+                 "no exact matching for '%s'. %d possible fields were found",
                  name, size(pos, 1));
         for j = 1:(size (pos, 1))
           fprintf ("%s\n", deblank (fields(pos(j),:)));
@@ -131,8 +128,7 @@
       endwhile
     endfor
 
-    if ((nargin == 2)
-        && (isstruct (varargin{2})))
+    if (nargin == 2 && isstruct (varargin{2}))
       ode_struct_value_check (varargin{2});
 
       optB_fields = fieldnames (varargin{2});
@@ -146,13 +142,14 @@
 
           if (size (pos, 1) == 0)
             warning ("OdePkg:InvalidArgument", ...
-                     "no property found with name ''%s''", name);
+                     "no property found with name '%s'", name);
           endif
 
           if (size(pos, 1) == 1)
-            if (! strcmp (lower (deblank (name)), lower (deblank (fields(pos,:)))))
+            if (! strcmp (lower (deblank (name)), ...
+                          lower (deblank (fields(pos,:)))))
               warning ("OdePkg:InvalidArgument", "no exact matching for ",
-                       "''%s''. Assuming you were intending ''%s''",
+                       "'%s'. Assuming you were intending '%s'",
                         name, deblank (fields(pos,:)));
             endif
             opt.(deblank (fields(pos,:))) = varargin{2}.(optB_fields{i});
@@ -160,7 +157,7 @@
           endif
 
           ## if there are more matching, ask the user to be more precise
-          warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ",
+          warning ("OdePkg:InvalidArgument", "no exact matching for '%s'. ",
                    "%d possible fields were found",
                    name, size (pos, 1));
           for j = 1:(size (pos, 1))
@@ -172,7 +169,7 @@
           until (ischar (name))
         endwhile
       endfor
-      return
+      return;
     endif
 
     ## if the second argument is not a struct,
@@ -198,12 +195,12 @@
 
         if (size (pos, 1) == 0)
           error ("OdePkg:InvalidArgument",
-                 "no property found with name ''%s''", name);
+                 "no property found with name '%s'", name);
         endif
 
         if (size (pos, 1) == 1)
           if (! strcmp (lower (deblank (name)), lower (deblank (fields(pos,:)))))
-            warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ",
+            warning ("OdePkg:InvalidArgument", "no exact matching for '%s'. ",
                      "%d possible fields were found",
                      name, size (pos, 1));
           endif
@@ -212,7 +209,7 @@
         endif
 
         ## if there are more matching, ask the user to be more precise
-        warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ",
+        warning ("OdePkg:InvalidArgument", "no exact matching for '%s'. ",
                  "%d possible fields were found",
                  name, size (pos, 1));
         for j = 1:(size (pos, 1))
@@ -227,7 +224,7 @@
 
     ## check if all has been done gives a valid OdePkg struct
     ode_struct_value_check (opt);
-    return
+    return;
   endif
 
   ## first input argument was not a struct
@@ -249,14 +246,14 @@
 
       if (size (pos, 1) == 0)
         error ("OdePkg:InvalidArgument",
-               "invalid property. No property found with name ''%s''", name);
+               "invalid property. No property found with name '%s'", name);
       endif
 
       if (size (pos, 1) == 1)
         if (! strcmp (lower (deblank (name)),
                       lower (deblank (fields(pos,:)))))
           warning ("OdePkg:InvalidArgument", "no exact matching for ",
-                   "''%s''. Assuming you were intending ''%s''",
+                   "'%s'. Assuming you were intending '%s'",
                    name, deblank (fields(pos,:)));
         endif
         opt.(deblank (fields(pos,:))) = varargin{i+1};
@@ -264,7 +261,7 @@
       endif
 
       ## if there are more matching, ask the user to be more precise
-      warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ",
+      warning ("OdePkg:InvalidArgument", "no exact matching for '%s'. ",
                "%d possible fields were found",
                name, size (pos, 1));
       for j = 1:(size (pos, 1))
@@ -279,38 +276,39 @@
 
   ## check if all has been done gives a valid OdePkg struct
   ode_struct_value_check (opt);
+
 endfunction
 
 ## function useful to print all the possible options
 function print_options ()
   
-  fprintf ("These following are all possible options.\n",
-           "Default values are put in square brackets.\n\n");
+  printf ("These following are all possible options.\n",
+          "Default values are put in square brackets.\n\n");
   
   disp ("             AbsTol:  scalar or vector, >0, [1.e-6]");
-  disp ("          Algorithm:  string, {[''gmres''], ''pcg'', ''bicgstab''}");
-  disp ("                BDF:  binary, {''on'', [''off'']}");
+  disp ("          Algorithm:  string, {['gmres'], 'pcg', 'bicgstab'}");
+  disp ("                BDF:  binary, {'on', ['off']}");
   disp ("             Choice:  switch, {[1], 2}");
   disp ("                Eta:  scalar, >=0, <1, [0.5]");
   disp ("             Events:  function_handle, []");
-  disp ("           Explicit:  binary, {''yes'', [''no'']}");
-  disp ("      InexactSolver:  string, {''inexact_newton'', ''fsolve'', []}");
+  disp ("           Explicit:  binary, {'yes', ['no']}");
+  disp ("      InexactSolver:  string, {'inexact_newton', 'fsolve', []}");
   disp ("       InitialSlope:  vector, []");
   disp ("        InitialStep:  scalar, >0, []");
   disp ("           Jacobian:  matrix or function_handle, []");
-  disp ("          JConstant:  binary, {''on'', [''off'']}");
+  disp ("          JConstant:  binary, {'on', ['off']}");
   disp ("           JPattern:  sparse matrix, []");
   disp ("               Mass:  matrix or function_handle, []");
-  disp ("       MassConstant:  binary, {''on'', [''off'']}");
-  disp ("       MassSingular:  switch, {''yes'', [''maybe''], ''no''}");
+  disp ("       MassConstant:  binary, {'on', ['off']}");
+  disp ("       MassSingular:  switch, {'yes', ['maybe'], 'no'}");
   disp ("MaxNewtonIterations:  scalar, integer, >0, [1.e3]");
   disp ("           MaxOrder:  switch, {1, 2, 3, 4, [5]}");
   disp ("            MaxStep:  scalar, >0, []");
-  disp ("   MStateDependence:  switch, {''none'', [''weak''], ''strong''}");
+  disp ("   MStateDependence:  switch, {'none', ['weak'], 'strong'}");
   disp ("          MvPattern:  sparse matrix, []");
   disp ("          NewtonTol:  scalar or vector, >0, []");
   disp ("        NonNegative:  vector of integers, []");
-  disp ("        NormControl:  binary, {''on'', [''off'']}");
+  disp ("        NormControl:  binary, {'on', ['off']}");
   disp ("          OutputFcn:  function_handle, []");
   disp ("         OutputSave:  scalar, integer, >0, []");
   disp ("          OutputSel:  scalar or vector, []");
@@ -319,30 +317,14 @@
   disp ("             Refine:  scalar, integer, >0, []");
   disp ("             RelTol:  scalar, >0, [1.e-3]");
   disp ("            Restart:  scalar, integer, >0, [20]");
-  disp ("              Stats:  binary, {''on'', [''off'']}");
+  disp ("              Stats:  binary, {'on', ['off']}");
   disp ("     TimeStepNumber:  scalar, integer, >0, []");
   disp ("       TimeStepSize:  scalar, >0, []");
-  disp ("        UseJacobian:  binary, {''yes'', [''no'']}");
-  disp ("         Vectorized:  binary, {''on'', [''off'']}");
+  disp ("        UseJacobian:  binary, {'yes', ['no']}");
+  disp ("         Vectorized:  binary, {'on', ['off']}");
 
 endfunction
 
-## All tests that are needed to check if a correct resp. valid option
-## has been set are implemented in ode_struct_value_check.m.
-%! ## Turn off output of warning messages for all tests, turn them on
-%! ## again if the last test is called
-%!  warning ('off', 'OdePkg:InvalidArgument');
-%!test odeoptA = odeset ();
-%!test odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1);
-%!     if (odeoptB.AbsTol != 1e-2), error; endif
-%!     if (odeoptB.RelTol != 1e-1), error; endif
-%!test odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1);
-%!     odeoptC = odeset (odeoptB, 'NormControl', 'on');
-%!test odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1);
-%!     odeoptC = odeset (odeoptB, 'NormControl', 'on');
-%!     odeoptD = odeset (odeoptC, odeoptB);
-%!
-%!  warning ('on', 'OdePkg:InvalidArgument');
 
 %!demo
 %! # A new OdePkg options structure with default values is created.
@@ -353,15 +335,29 @@
 %! # A new OdePkg options structure with manually set options 
 %! # "AbsTol" and "RelTol" is created.
 %!
-%! odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1);
+%! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
 %!
 %!demo
 %! # A new OdePkg options structure from odeoptB is created with
 %! # a modified value for option "NormControl".
 %!
-%! odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1);
-%! odeoptC = odeset (odeoptB, 'NormControl', 'on');
+%! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
+%! odeoptC = odeset (odeoptB, "NormControl", "on");
 
-## Local Variables: ***
-## mode: octave ***
-## End: ***
+## All tests that are needed to check if a correct resp. valid option
+## has been set are implemented in ode_struct_value_check.m.
+%! ## Turn off output of warning messages for all tests, turn them on
+%! ## again if the last test is called
+%!  warning ("off", "OdePkg:InvalidArgument");
+%!test odeoptA = odeset ();
+%!test odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
+%!     if (odeoptB.AbsTol != 1e-2), error; endif
+%!     if (odeoptB.RelTol != 1e-1), error; endif
+%!test odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
+%!     odeoptC = odeset (odeoptB, "NormControl", "on");
+%!test odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
+%!     odeoptC = odeset (odeoptB, "NormControl", "on");
+%!     odeoptD = odeset (odeoptC, odeoptB);
+%!
+%!  warning ("on", "OdePkg:InvalidArgument");
+