--- a/scripts/ode/ode23.m	Wed Oct 05 23:09:31 2016 -0400
+++ b/scripts/ode/ode23.m	Thu Oct 06 06:54:32 2016 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
 ## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com>
 ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
 ## Copyright (C) 2006-2016 Thomas Treichl <treichl@users.sourceforge.net>
@@ -45,10 +46,7 @@
 ## By default, @code{ode23} 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.
+## default value of @math{1e-6}. 
 ##
 ## @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
@@ -106,25 +104,25 @@
     print_usage ();
   endif
 
-  order = 3;
+  order  = 3;
   solver = "ode23";
 
   if (nargin >= 4)
     if (! isstruct (varargin{1}))
       ## varargin{1:len} are parameters for fun
       odeopts = odeset ();
-      odeopts.funarguments = varargin;
+      funarguments = varargin;
     elseif (length (varargin) > 1)
       ## varargin{1} is an ODE options structure opt
-      odeopts = ode_struct_value_check ("ode23", varargin{1}, "ode23");
-      odeopts.funarguments = {varargin{2:length(varargin)}};
+      odeopts = varargin{1};
+      funarguments = {varargin{2:length(varargin)}};
     else  # if (isstruct (varargin{1}))
-      odeopts = ode_struct_value_check ("ode23", varargin{1}, "ode23");
-      odeopts.funarguments = {};
+      odeopts = varargin{1};
+      funarguments = {};
     endif
   else  # nargin == 3
     odeopts = odeset ();
-    odeopts.funarguments = {};
+    funarguments = {};
   endif
 
   if (! isnumeric (trange) || ! isvector (trange))
@@ -132,17 +130,14 @@
            "ode23: TRANGE must be a numeric vector");
   endif
 
-  TimeStepNumber = odeget (odeopts, "TimeStepNumber", [], "fast");
-  TimeStepSize = odeget (odeopts, "TimeStepSize", [], "fast");
-  if (length (trange) < 2
-      && (isempty (TimeStepSize) || isempty (TimeStepNumber)))
+  if (length (trange) < 2)
     error ("Octave:invalid-input-arg",
            "ode23: TRANGE must contain at least 2 elements");
   elseif (trange(2) == trange(1))
     error ("Octave:invalid-input-arg",
            "ode23: invalid time span, TRANGE(1) == TRANGE(2)");
   else
-    odeopts.direction = sign (trange(2) - trange(1));
+    direction = sign (trange(2) - trange(1));
   endif
   trange = trange(:);
 
@@ -164,54 +159,33 @@
            "ode23: FUN must be a valid function handle");
   endif
 
+
   ## 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";
-    odeopts.stepsizefixed = false;
-  elseif (! isempty (TimeStepNumber) && ! isempty (TimeStepSize))
-    integrate_func = "n_steps";
-    odeopts.stepsizefixed = true;
-    if (sign (TimeStepSize) != odeopts.direction)
-      warning ("Octave:invalid-input-arg",
-               ["ode23: option \"TimeStepSize\" has the wrong sign, ", ...
-                "but will be corrected automatically\n"]);
-      TimeStepSize = -TimeStepSize;
-    endif
-  elseif (isempty (TimeStepNumber) && ! isempty (TimeStepSize))
-    integrate_func = "const";
-    odeopts.stepsizefixed = true;
-    if (sign (TimeStepSize) != odeopts.direction)
-      warning ("Octave:invalid-input-arg",
-               ["ode23: option \"TimeStepSize\" has the wrong sign, ",
-                "but will be corrected automatically\n"]);
-      TimeStepSize = -TimeStepSize;
-    endif
-  else
-    warning ("Octave:invalid-input-arg",
-             "ode23: assuming an adaptive integrate function\n");
-    integrate_func = "adaptive";
-  endif
+
+
+  persistent defaults   = [];
+  persistent classes    = [];
+  persistent attributes = [];
+
+
+  [defaults, classes, attributes] = odedefaults (numel (init), trange(1),
+                                                 trange(end));
 
-  if (isempty (odeopts.RelTol) && ! odeopts.stepsizefixed)
-    odeopts.RelTol = 1e-3;
-  elseif (! isempty (odeopts.RelTol) && odeopts.stepsizefixed)
-    warning ("Octave:invalid-input-arg",
-             ["ode23: option \"RelTol\" is ignored", ...
-              " when fixed time stamps are given\n"]);
-  endif
+  defaults   = rmfield (defaults,   {"Jacobian", "JPattern", "Vectorized", ...
+                                     "MvPattern", "MassSingular", ...
+                                     "InitialSlope", "MaxOrder", "BDF"});
+  classes    = rmfield (classes,    {"Jacobian", "JPattern", "Vectorized", ...
+                                     "MvPattern", "MassSingular", ...
+                                     "InitialSlope", "MaxOrder", "BDF"});
+  attributes = rmfield (attributes, {"Jacobian", "JPattern", "Vectorized", ...
+                                     "MvPattern", "MassSingular", ...
+                                     "InitialSlope", "MaxOrder", "BDF"});
 
-  if (isempty (odeopts.AbsTol) && ! odeopts.stepsizefixed)
-    odeopts.AbsTol = 1e-6;
-  elseif (! isempty (odeopts.AbsTol) && odeopts.stepsizefixed)
-    warning ("Octave:invalid-input-arg",
-             ["ode23: option \"AbsTol\" is ignored", ...
-              " when fixed time stamps are given\n"]);
-  else
-    odeopts.AbsTol = odeopts.AbsTol(:);  # Create column vector
-  endif
+  odeopts = odemergeopts (odeopts, defaults, classes, attributes, 'ode23');
 
-  odeopts.normcontrol = strcmp (odeopts.NormControl, "on");
+  odeopts.funarguments = funarguments;
+  odeopts.direction    = direction;
 
   if (! isempty (odeopts.NonNegative))
     if (isempty (odeopts.Mass))
@@ -233,47 +207,15 @@
     odeopts.haveoutputfunction = ! isempty (odeopts.OutputFcn);
   endif
 
-  odeopts.haveoutputselection = ! isempty (odeopts.OutputSel);
-
-  if (odeopts.Refine > 0)
-    odeopts.haverefine = true;
-  else
-    odeopts.haverefine = false;
-  endif
-
-  if (isempty (odeopts.InitialStep) && strcmp (integrate_func, "adaptive"))
+  if (isempty (odeopts.InitialStep))
     odeopts.InitialStep = odeopts.direction * ...
                           starting_stepsize (order, fun, trange(1),
                                              init, odeopts.AbsTol,
                                              odeopts.RelTol,
-                                             odeopts.normcontrol);
-  elseif (isempty (odeopts.InitialStep))
-    odeopts.InitialStep = TimeStepSize;
-  endif
-
-  if (isempty (odeopts.MaxStep) && ! odeopts.stepsizefixed)
-    odeopts.MaxStep = abs (trange(end) - trange(1)) / 10;
+                                             strcmp (odeopts.NormControl,
+                                             "on"));
   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 (! isempty (odeopts.Jacobian))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"Jacobian\" is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.JPattern))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"JPattern\" is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.Vectorized))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"Vectorized\" is ignored by this solver\n");
-  endif
 
   if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass))
     havemasshandle = false;
@@ -284,75 +226,38 @@
     havemasshandle = false; # mass = diag (ones (length (init), 1), 0);
   endif
 
-  massdependence = ! strcmp (odeopts.MStateDependence, "none");
-
-  ## Other options that are not used by this solver.
-  if (! isempty (odeopts.MvPattern))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"MvPattern\" is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.MassSingular))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"MassSingular\" is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.InitialSlope))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"InitialSlope\" is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.MaxOrder))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"MaxOrder\" is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.BDF))
-    warning ("Octave:invalid-input-arg",
-             "ode23: option \"BDF\" is ignored by this solver\n");
-  endif
 
   ## Starting the initialization of the core solver ode23
 
   if (havemasshandle)   # Handle only the dynamic mass matrix,
-    if (massdependence) # constant mass matrices have already
+    if (! strcmp (odeopts.MStateDependence, "none")) # 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)
+    else                 # if ((! strcmp (odeopts.MStateDependence, "none")) == 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_23, ...
-                                     order, fun, trange, init, odeopts);
-    case "n_steps"
-      solution = integrate_n_steps (@runge_kutta_23, ...
-                                    fun, trange(1), init, ...
-                                    TimeStepSize, TimeStepNumber, odeopts);
-    case "const"
-      solution = integrate_const (@runge_kutta_23, ...
-                                  fun, trange, init, ...
-                                  TimeStepSize, odeopts);
-  endswitch
+  
+  solution = integrate_adaptive (@runge_kutta_23, ...
+                                 order, fun, trange, init, odeopts);
+    
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
     feval (odeopts.OutputFcn, solution.t(end), ...
            solution.x(end,:)', "done", odeopts.funarguments{:});
   endif
-  if (odeopts.haveeventfunction)   # Cleanup event function handling
+  if (! isempty (odeopts.Events))   # Cleanup event function handling
     ode_event_handler (odeopts.Events, solution.t(end), ...
-                         solution.x(end,:)', "done", odeopts.funarguments{:});
+                       solution.x(end,:)', "done", odeopts.funarguments{:});
   endif
 
   ## Print additional information if option Stats is set
   if (strcmp (odeopts.Stats, "on"))
-    havestats = true;
     nsteps    = solution.cntloop;             # cntloop from 2..end
     nfailed   = solution.cntcycles - nsteps;  # cntcycl from 1..end
     nfevals   = 3 * solution.cntcycles + 1;   # number of ode evaluations
@@ -365,8 +270,6 @@
       printf ("Number of failed attempts:  %d\n", nfailed);
       printf ("Number of function calls:   %d\n", nfevals);
     endif
-  else
-    havestats = false;
   endif
 
   if (nargout == 2)
@@ -376,12 +279,12 @@
     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 = solver;   # Solver name is saved in field solver
-    if (odeopts.haveeventfunction)
+    if (! isempty (odeopts.Events))
       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 (havestats)
+    if (strcmp (odeopts.Stats, "on"))
       varargout{1}.stats = struct ();
       varargout{1}.stats.nsteps   = nsteps;
       varargout{1}.stats.nfailed  = nfailed;
@@ -394,7 +297,7 @@
     varargout = cell (1,5);
     varargout{1} = solution.t;
     varargout{2} = solution.x;
-    if (odeopts.haveeventfunction)
+    if (! isempty (odeopts.Events))
       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
@@ -489,10 +392,6 @@
 %! opt = odeset;
 %! [t, y] = ode23 (@fpol, [0 2], [2 0], opt, 12, 13, "KL");
 %! assert ([t(end), y(end,:)], [2, fref], 1e-2);
-%!test  # Solve vdp in fixed step sizes
-%! opt = odeset("TimeStepSize",0.1);
-%! [t, y] = ode23 (@fpol, [0,2], [2 0], opt);
-%! assert (t(:), [0:0.1:2]', 1e-3);
 %!test  # Solve another anonymous function below zero
 %! ref = [0, 14.77810590694212];
 %! [t, y] = ode23 (@(t,y) y, [-2 0], 2);
@@ -612,7 +511,4 @@
 %!  ode23 (@fpol, [0 25], [3 15 1; 3 15 1]);
 %!error <FUN must be a valid function handle>
 %!  ode23 (1, [0 25], [3 15 1]);
-%!error  # strange ODEOPT structure
-%!  opt = struct ("foo", 1);
-%!  [t, y] = ode23 (@fpol, [0 2], [2 0], opt);
 

--- a/scripts/ode/ode45.m	Wed Oct 05 23:09:31 2016 -0400
+++ b/scripts/ode/ode45.m	Thu Oct 06 06:54:32 2016 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
 ## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com>
 ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
 ## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net>
@@ -43,10 +44,7 @@
 ## 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.
+## default value of @math{1e-6}.
 ##
 ## @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
@@ -89,25 +87,25 @@
     print_usage ();
   endif
 
-  order = 5;  # runge_kutta_45_dorpri uses local extrapolation
+  order  = 5;  # runge_kutta_45_dorpri uses local extrapolation
   solver = "ode45";
 
   if (nargin >= 4)
     if (! isstruct (varargin{1}))
       ## varargin{1:len} are parameters for fun
       odeopts = odeset ();
-      odeopts.funarguments = varargin;
+      funarguments = varargin;
     elseif (length (varargin) > 1)
-      ## varargin{1} is an ODE options structure vopt
-      odeopts = ode_struct_value_check ("ode45", varargin{1}, "ode45");
-      odeopts.funarguments = {varargin{2:length(varargin)}};
+      ## varargin{1} is an ODE options structure opt
+      odeopts = varargin{1};
+      funarguments = {varargin{2:length(varargin)}};
     else  # if (isstruct (varargin{1}))
-      odeopts = ode_struct_value_check ("ode45", varargin{1}, "ode45");
-      odeopts.funarguments = {};
+      odeopts = varargin{1};
+      funarguments = {};
     endif
   else  # nargin == 3
     odeopts = odeset ();
-    odeopts.funarguments = {};
+    funarguments = {};
   endif
 
   if (! isnumeric (trange) || ! isvector (trange))
@@ -115,17 +113,15 @@
            "ode45: TRANGE must be a numeric vector");
   endif
 
-  TimeStepNumber = odeget (odeopts, "TimeStepNumber", [], "fast");
-  TimeStepSize = odeget (odeopts, "TimeStepSize", [], "fast");
-  if (length (trange) < 2
-      && (isempty (TimeStepSize) || isempty (TimeStepNumber)))
+
+  if (length (trange) < 2)
     error ("Octave:invalid-input-arg",
            "ode45: TRANGE must contain at least 2 elements");
   elseif (trange(1) == trange(2))
     error ("Octave:invalid-input-arg",
            "ode45: invalid time span, TRANGE(1) == TRANGE(2)");
   else
-    odeopts.direction = sign (trange(2) - trange(1));
+    direction = sign (trange(2) - trange(1));
   endif
   trange = trange(:);
 
@@ -147,54 +143,33 @@
            "ode45: FUN must be a valid function handle");
   endif
 
+
   ## 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";
-    odeopts.stepsizefixed = false;
-  elseif (! isempty (TimeStepNumber) && ! isempty (TimeStepSize))
-    integrate_func = "n_steps";
-    odeopts.stepsizefixed = true;
-    if (sign (TimeStepSize) != odeopts.direction)
-      warning ("Octave:invalid-input-arg",
-               ["ode45: option 'TimeStepSize' has the wrong sign, ", ...
-                "but will be corrected automatically\n"]);
-      TimeStepSize = -TimeStepSize;
-    endif
-  elseif (isempty (TimeStepNumber) && ! isempty (TimeStepSize))
-    integrate_func = "const";
-    odeopts.stepsizefixed = true;
-    if (sign (TimeStepSize) != odeopts.direction)
-      warning ("Octave:invalid-input-arg",
-               ["ode45: option 'TimeStepSize' has the wrong sign, ",
-                "but will be corrected automatically\n"]);
-      TimeStepSize = -TimeStepSize;
-    endif
-  else
-    warning ("Octave:invalid-input-arg",
-             "ode45: assuming an adaptive integrate function\n");
-    integrate_func = "adaptive";
-  endif
+
+  persistent defaults   = [];
+  persistent classes    = [];
+  persistent attributes = [];
+
+
+  [defaults, classes, attributes] = odedefaults (numel (init), trange(1),
+                                                 trange(end));
 
-  if (isempty (odeopts.RelTol) && ! odeopts.stepsizefixed)
-    odeopts.RelTol = 1e-3;
-  elseif (! isempty (odeopts.RelTol) && odeopts.stepsizefixed)
-    warning ("Octave:invalid-input-arg",
-             ["ode45: option 'RelTol' is ignored", ...
-              " when fixed time stamps are given\n"]);
-  endif
+  defaults   = odeset (defaults, 'Refine', 4);
+  defaults   = rmfield (defaults,   {"Jacobian", "JPattern", "Vectorized", ...
+                                     "MvPattern", "MassSingular", ...
+                                     "InitialSlope", "MaxOrder", "BDF"});
+  classes    = rmfield (classes,    {"Jacobian", "JPattern", "Vectorized", ...
+                                     "MvPattern", "MassSingular", ...
+                                     "InitialSlope", "MaxOrder", "BDF"});
+  attributes = rmfield (attributes, {"Jacobian", "JPattern", "Vectorized", ...
+                                     "MvPattern", "MassSingular", ...
+                                     "InitialSlope", "MaxOrder", "BDF"}); 
 
-  if (isempty (odeopts.AbsTol) && ! odeopts.stepsizefixed)
-    odeopts.AbsTol = 1e-6;
-  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
+  odeopts = odemergeopts (odeopts, defaults, classes, attributes, 'ode45');
 
-  odeopts.normcontrol = strcmp (odeopts.NormControl, "on");
+  odeopts.funarguments = funarguments;
+  odeopts.direction    = direction;
 
   if (! isempty (odeopts.NonNegative))
     if (isempty (odeopts.Mass))
@@ -216,48 +191,15 @@
     odeopts.haveoutputfunction = ! isempty (odeopts.OutputFcn);
   endif
 
-  odeopts.haveoutputselection = ! isempty (odeopts.OutputSel);
-
-  if (odeopts.Refine > 0)
-    odeopts.haverefine = true;
-  else
-    odeopts.haverefine = false;
-  endif
-
-  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);
-  elseif (isempty (odeopts.InitialStep))
-    odeopts.InitialStep = TimeStepSize;
-  endif
+  if (isempty (odeopts.InitialStep))
+    odeopts.InitialStep = odeopts.direction * ...
+                          starting_stepsize (order, fun, trange(1),
+                                             init, odeopts.AbsTol,
+                                             odeopts.RelTol,
+                                             strcmp (odeopts.NormControl,
+                                             "on"));
+  endif 
 
-  if (isempty (odeopts.MaxStep) && ! odeopts.stepsizefixed)
-    odeopts.MaxStep = abs (trange(end) - trange(1)) / 10;
-  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 (! isempty (odeopts.Jacobian))
-    warning ("Octave:invalid-input-arg",
-             "ode45: option 'Jacobian' is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.JPattern))
-    warning ("Octave:invalid-input-arg",
-             "ode45: option 'JPattern' is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.Vectorized))
-    warning ("Octave:invalid-input-arg",
-             "ode45: option 'Vectorized' is ignored by this solver\n");
-  endif
 
   if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass))
     havemasshandle = false;
@@ -268,76 +210,38 @@
     havemasshandle = false;   # mass = diag (ones (length (init), 1), 0);
   endif
 
-  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
-
-  if (! isempty (odeopts.MassSingular))
-    warning ("Octave:invalid-input-arg",
-             "ode45: option 'MassSingular' is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.InitialSlope))
-    warning ("Octave:invalid-input-arg",
-             "ode45: option 'InitialSlope' is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.MaxOrder))
-    warning ("Octave:invalid-input-arg",
-             "ode45: option 'MaxOrder' is ignored by this solver\n");
-  endif
-
-  if (! isempty (odeopts.BDF))
-    warning ("Octave:invalid-input-arg",
-             "ode45: option 'BDF' is ignored by this solver\n");
-  endif
 
   ## Starting the initialization of the core solver ode45
 
   if (havemasshandle)   # Handle only the dynamic mass matrix,
-    if (massdependence) # constant mass matrices have already
+    if (! strcmp (odeopts.MStateDependence, "none")) # 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)
+    else                 # if ((! strcmp (odeopts.MStateDependence, "none")) == 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,
-                                     order, fun, trange, init, odeopts);
-    case "n_steps"
-      solution = integrate_n_steps (@runge_kutta_45_dorpri,
-                                    fun, trange(1), init,
-                                    TimeStepSize, TimeStepNumber, odeopts);
-    case "const"
-      solution = integrate_const (@runge_kutta_45_dorpri,
-                                  fun, trange, init,
-                                  TimeStepSize, odeopts);
-  endswitch
+ 
+  solution = integrate_adaptive (@runge_kutta_45_dorpri,
+                                 order, fun, trange, init, odeopts);
+  
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, solution.t(end),
+    feval (odeopts.OutputFcn, solution.t(end), ...
            solution.x(end,:)', "done", odeopts.funarguments{:});
   endif
-  if (odeopts.haveeventfunction)   # Cleanup event function handling
-    ode_event_handler (odeopts.Events, solution.t(end),
-                         solution.x(end,:)', "done",
-                         odeopts.funarguments{:});
+  if (! isempty (odeopts.Events))   # Cleanup event function handling
+    ode_event_handler (odeopts.Events, solution.t(end), ...
+                       solution.x(end,:)', "done", odeopts.funarguments{:});
   endif
 
   ## Print additional information if option Stats is set
   if (strcmp (odeopts.Stats, "on"))
-    havestats = true;
     nsteps    = solution.cntloop;             # cntloop from 2..end
     nfailed   = solution.cntcycles - nsteps;  # cntcycl from 1..end
     nfevals   = 6 * solution.cntcycles + 1;   # number of ode evaluations
@@ -350,8 +254,6 @@
       printf ("Number of failed attempts:  %d\n", nfailed);
       printf ("Number of function calls:   %d\n", nfevals);
     endif
-  else
-    havestats = false;
   endif
 
   if (nargout == 2)
@@ -361,12 +263,12 @@
     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 = solver;   # Solver name is saved in field solver
-    if (odeopts.haveeventfunction)
+    if (! isempty (odeopts.Events))
       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 (havestats)
+    if (strcmp (odeopts.Stats, "on"))
       varargout{1}.stats = struct ();
       varargout{1}.stats.nsteps   = nsteps;
       varargout{1}.stats.nfailed  = nfailed;
@@ -379,7 +281,7 @@
     varargout = cell (1,5);
     varargout{1} = solution.t;
     varargout{2} = solution.x;
-    if (odeopts.haveeventfunction)
+    if (! isempty (odeopts.Events))
       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
@@ -480,10 +382,6 @@
 %! opt = odeset;
 %! [t, y] = ode45 (@fpol, [0 2], [2 0], opt, 12, 13, "KL");
 %! assert ([t(end), y(end,:)], [2, fref], 1e-2);
-%!test  # Solve vdp in fixed step sizes
-%! 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];
 %! [t, y] = ode45 (@(t,y) y, [-2 0], 2);
@@ -609,7 +507,4 @@
 %!  ode45 (@fpol, [0 25], [3 15 1; 3 15 1]);
 %!error <FUN must be a valid function handle>
 %!  ode45 (1, [0 25], [3 15 1]);
-%!error  # strange ODEOPT structure
-%! opt = struct ("foo", 1);
-%! [t, y] = ode45 (@fpol, [0 2], [2 0], opt);
 

--- a/scripts/ode/odeget.m	Wed Oct 05 23:09:31 2016 -0400
+++ b/scripts/ode/odeget.m	Thu Oct 06 06:54:32 2016 +0200
@@ -1,3 +1,5 @@
+## Copyright (C) 2016, Carlo de Falco
+## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
 ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
 ## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net>
 ##
@@ -39,79 +41,19 @@
 
 function val = odeget (ode_opt, field, default = [], opt = "")
 
-  ## Shortcut for quickly extracting option
-  if (strncmp (opt, "fast", 4))
-    try
-      val = ode_opt.(field);
-      if (strcmp (opt, "fast_not_empty") && isempty (val))
-        val = default;
-      endif
-    catch
-      val = default;
-    end_try_catch
-    return;
-  endif
-
-  if (nargin < 1 || nargin > 4)
-    print_usage ();
-  endif
-
-  ## Shortcut for empty option structures
-  if (isempty (ode_opt))
-    if (nargin < 3)
-      val = [];
-    else
+  validateattributes (ode_opt, {'struct'}, {'nonempty'});
+  validateattributes (field, {'char'}, {'nonempty'});
+  
+  if (! isfield (ode_opt, field))
+    error ('Octave:odeget:InvalidPropName',
+           'odeget: Unrecognized property name "%s".', field);
+  else
+    val = ode_opt.(field);
+    if (isempty (val)) 
       val = default;
     endif
-    return;
   endif
-
-  if (! isstruct (ode_opt))
-    error ("odeget: ODE_OPT must be a valid ODE_STRUCT");
-  elseif (! ischar (field))
-    error ("odeget: FIELD must be a string");
-  endif
-
-  ## Check if the given struct is a valid ODEOPT struct
-  ode_struct_value_check ("odeget", ode_opt);
-
-  ## Define all the possible ODEOPT fields
-  persistent options = known_option_names ();
-
-  exactmatch = true;
-  match = find (strcmpi (field, options));
-  if (isempty (match))
-    match = find (strncmpi (field, options, length (field)));
-    exactmatch = false;
-  endif
-
-  if (isempty (match))
-    ## Possibly a custom user-defined option
-    try
-      val = ode_opt.(field);
-    catch
-      warning ("Octave:invalid-input-arg",
-               "odeget: no field '%s' exists in ODE_OPT\n", field);
-      val = default;
-    end_try_catch
-  elseif (numel (match) == 1)
-    if (! exactmatch)
-      warning ("odeget:NoExactMatching",
-               "odeget: no exact match for '%s'.  Assuming '%s'.\n",
-               field, options{match});
-    endif
-    val = [];
-    try
-      val = ode_opt.(options{match});
-    end_try_catch
-    if (isempty (val))
-      val = default;
-    endif
-  else
-    error ("odeget: no exact match for '%s'.  Possible fields found: %s.",
-           field, strjoin (options(match), ", "));
-  endif
-
+  
 endfunction
 
 
@@ -129,16 +71,14 @@
 %!assert (odeget (odeset (), "Stats"), [])
 %!assert (odeget (odeset (), "Stats", "on"), "on")
 %!assert (odeget (odeset (), "Mass"), [])
-%!assert (odeget (odeset (), "AbsTol", 1e-6, "fast"), [])
-%!assert (odeget (odeset (), "AbsTol", 1e-6, "fast_not_empty"), 1e-6)
 %!assert (odeget (odeset (), "AbsTol", 1e-9), 1e-9)
+%!assert (odeget (odeset ("AbsTol", 1e-9), "AbsTol", []), 1e-9)
+%!assert (odeget (odeset ('foo', 42), 'foo'), 42)
 
 %!error odeget ()
 %!error odeget (1)
 %!error odeget (1,2,3,4,5)
-%!error <ODE_OPT must be a valid ODE_STRUCT> odeget (1, "opt1")
-%!error <FIELD must be a string> odeget (struct ("opt1", 1), 1)
-%!warning <no field 'foo' exists> odeget (struct ("opt1", 1), "foo");
-%!warning <no exact match for 'Rel'.  Assuming 'RelTol'> odeget (struct ("RelTol", 1), "Rel");
-%!error <Possible fields found: InitialSlope, InitialStep> odeget (odeset (), "Initial")
+%!error odeget (1, "opt1")
+%!error odeget (struct ("opt1", 1), 1)
+%!error odeget (struct ("opt1", 1), "foo");
 

--- a/scripts/ode/odeset.m	Wed Oct 05 23:09:31 2016 -0400
+++ b/scripts/ode/odeset.m	Thu Oct 06 06:54:32 2016 +0200
@@ -1,3 +1,5 @@
+## Copyright (C) 2016, Carlo de Falco
+## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
 ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
 ## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net>
 ##
@@ -47,168 +49,59 @@
 
 function odestruct = odeset (varargin)
 
-  ## Column vector of all possible ODE options
-  persistent options = known_option_names ();
+  persistent p;
+
+  if (isempty (p))
 
-  if (nargin == 0)
-    ## Special calling syntax to display defaults
-    if (nargout == 0)
-      print_options ();
-    else
-      odestruct = cell2struct (cell (numel (options), 1), options);
-    endif
-    return;
+    p = inputParser ();
+    p.addParameter ("AbsTol", []);
+    p.addParameter ("BDF", []);
+    p.addParameter ("Events", []);
+    p.addParameter ("InitialSlope", []);
+    p.addParameter ("InitialStep", []);
+    p.addParameter ("Jacobian", []);
+    p.addParameter ("JConstant", []);
+    p.addParameter ("JPattern", []);
+    p.addParameter ("Mass", []);
+    p.addParameter ("MassConstant", []);
+    p.addParameter ("MassSingular", []);
+    p.addParameter ("MaxOrder", []);
+    p.addParameter ("MaxStep", []);
+    p.addParameter ("MStateDependence", []);
+    p.addParameter ("MvPattern", []);
+    p.addParameter ("NonNegative", []);
+    p.addParameter ("NormControl", []);
+    p.addParameter ("OutputFcn", []);
+    p.addParameter ("OutputSel", []);
+    p.addParameter ("Refine", []);
+    p.addParameter ("RelTol", []);
+    p.addParameter ("Stats", []);
+    p.addParameter ("Vectorized", []);
+    p.KeepUnmatched = true;
+    
   endif
 
-  ## initialize output
-  odestruct = cell2struct (cell (numel (options), 1), options);
-
-  if (isstruct (varargin{1}))
-    oldstruct = varargin{1};
-
-    ## Copy oldstruct values into output odestruct
-    for [val, name] = oldstruct
-
-      exactmatch = true;
-      match = find (strcmpi (name, options));
-      if (isempty (match))
-        match = find (strncmpi (name, options, length (name)));
-        exactmatch = false;
-      endif
-
-      if (isempty (match))
-        odestruct.(name) = val;
-      elseif (numel (match) == 1)
-        if (! exactmatch)
-          warning ("odeset:NoExactMatching",
-                   "no exact match for '%s'.  Assuming '%s'.",
-                   name, options{match});
-        endif
-        odestruct.(options{match}) = val;
-      else
-        error ("odeset: no exact match for '%s'.  Possible fields found: %s.",
-               name, strjoin (options(match), ", "));
-      endif
-
-      if (nargin == 1)
-        ## Check if all changes have resulted in a valid ODEOPT struct
-        ode_struct_value_check ("odeset", odestruct);
-        return;
-      endif
+  if (nargin == 0 && nargout == 0)
+    print_options ();
+  else
+    p.parse (varargin{:});
+    odestruct = p.Results;
+    odestruct_extra = p.Unmatched;
 
-    endfor
-
-    ## At this point, odestruct has been initialized with default values,
-    ## and if oldstruct was present it has overwritten fields in odestruct.
-
-    if (nargin == 2 && isstruct (varargin{2}))
-      newstruct = varargin{2};
-
-      ## Update the first struct with the values from the second one
-      for [val, name] = newstruct
-
-        exactmatch = true;
-        match = find (strcmpi (name, options));
-        if (isempty (match))
-          match = find (strncmpi (name, options, length (name)));
-          exactmatch = false;
-        endif
-
-        if (isempty (match))
-          odestruct.(name) = val;
-        elseif (numel (match) == 1)
-          if (! exactmatch)
-            warning ("odeset:NoExactMatching",
-                     "no exact match for '%s'.  Assuming '%s'.",
-                     name, options{match});
-          endif
-          odestruct.(options{match}) = val;
-        else
-          error ("odeset: no exact match for '%s'.  Possible fields found: %s.",
-                 name, strjoin (options(match), ", "));
-        endif
-      endfor
-
-      ## Check if all changes have resulted in a valid ODEOPT struct
-      ode_struct_value_check ("odeset", odestruct);
-      return;
-    endif
+    s1 = cellfun (@(x) ifelse (iscell(x), {x}, x),
+                  struct2cell(odestruct),
+                  'UniformOutput', false);
 
-    ## Second argument is not a struct
-    if (mod (nargin, 2) != 1)
-      error ("odeset: FIELD/VALUE arguments must occur in pairs");
-    endif
-    if (! all (cellfun ("isclass", varargin(2:2:end), "char")))
-      error ("odeset: All FIELD names must be strings");
-    endif
-
-    ## Write new field/value pairs into odestruct
-    for i = 2:2:nargin
-      name = varargin{i};
-
-      exactmatch = true;
-      match = find (strcmpi (name, options));
-      if (isempty (match))
-        match = find (strncmpi (name, options, length (name)));
-        exactmatch = false;
-      endif
-
-      if (isempty (match))
-        odestruct.(name) = varargin{i+1};
-      elseif (numel (match) == 1)
-        if (! exactmatch)
-          warning ("odeset:NoExactMatching",
-                   "no exact match for '%s'.  Assuming '%s'.",
-                   name, options{match});
-        endif
-        odestruct.(options{match}) = varargin{i+1};
-      else
-        error ("odeset: no exact match for '%s'.  Possible fields found: %s.",
-               name, strjoin (options(match), ", "));
-      endif
-    endfor
-
-    ## Check if all changes have resulted in a valid ODEOPT struct
-    ode_struct_value_check ("odeset", odestruct);
-
-  else
-    ## First input argument was not a struct, must be field/value pairs
-    if (mod (nargin, 2) != 0)
-      error ("odeset: FIELD/VALUE arguments must occur in pairs");
-    elseif (! all (cellfun ("isclass", varargin(1:2:end), "char")))
-      error ("odeset: All FIELD names must be strings");
-    endif
-
-    for i = 1:2:nargin
-      name = varargin{i};
-
-      exactmatch = true;
-      match = find (strcmpi (name, options));
-      if (isempty (match))
-        match = find (strncmpi (name, options, length (name)));
-        exactmatch = false;
-      endif
-
-      if (isempty (match))
-        odestruct.(name) = varargin{i+1};
-      elseif (numel (match) == 1)
-        if (! exactmatch)
-          warning ("odeset:NoExactMatching",
-                   "no exact match for '%s'.  Assuming '%s'.",
-                   name, options{match});
-        endif
-        odestruct.(options{match}) = varargin{i+1};
-      else
-        error ("odeset: no exact match for '%s'.  Possible fields found: %s.",
-               name, strjoin (options(match), ", "));
-      endif
-    endfor
-
-    ## Check if all changes have resulted in a valid ODEOPT struct
-    ode_struct_value_check ("odeset", odestruct);
-
+    s2 = cellfun (@(x) ifelse (iscell(x), {x}, x),
+                  struct2cell(odestruct_extra),
+                  'UniformOutput', false);
+    
+    C = [fieldnames(odestruct)       s1;
+         fieldnames(odestruct_extra) s2];
+    
+    odestruct = struct (C'{:});
   endif
-
+  
 endfunction
 
 ## function to print all possible options
@@ -262,11 +155,7 @@
 %! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
 %! odeoptC = odeset (odeoptB, "NormControl", "on");
 
-## All tests that are needed to check if a valid option has been set are
-## implemented in ode_struct_value_check.m
-## FIXME: xtest currently fails as there are two extra options to control
-##        fixed step integration options.
-%!xtest
+%!test
 %! odeoptA = odeset ();
 %! assert (isstruct (odeoptA));
 %! assert (numfields (odeoptA), 23);
@@ -298,10 +187,13 @@
 %! end_unwind_protect
 
 ## Test input validation
-%!error <FIELD/VALUE arguments must occur in pairs> odeset ("opt1")
-%!error <FIELD names must be strings> odeset (1, 1)
-%!error <FIELD/VALUE arguments must occur in pairs> odeset (odeset (), "opt1")
-%!error <FIELD names must be strings> odeset (odeset (), 1, 1)
-%!warning <no exact match for 'Rel'.  Assuming 'RelTol'> odeset ("Rel", 1);
-%!error <Possible fields found: InitialSlope, InitialStep> odeset ("Initial", 1)
+%!error <argument 'OPT1' is not a valid parameter> odeset ("opt1")
+%!error  odeset (1, 1)
+%!error <argument 'OPT1' is not a valid parameter> odeset (odeset (), "opt1")
+%!error  odeset (odeset (), 1, 1)
 
+##FIXME: Add not exact match option 
+## %!warning <no exact match for 'Rel'.  Assuming 'RelTol'> odeset ("Rel", 1);
+## %!error <Possible fields found: InitialSlope, InitialStep> odeset ("Initial", 1)
+
+

--- a/scripts/ode/private/integrate_adaptive.m	Wed Oct 05 23:09:31 2016 -0400
+++ b/scripts/ode/private/integrate_adaptive.m	Thu Oct 06 06:54:32 2016 +0200
@@ -76,7 +76,7 @@
   dt = odeget (options, "InitialStep", [], "fast");
   if (isempty (dt))
     dt = starting_stepsize (order, func, t, x, options.AbsTol, options.RelTol,
-                            options.normcontrol);
+                            strcmp (options.NormControl, "on"));
   endif
 
   dir = odeget (options, "direction", [], "fast");
@@ -91,7 +91,7 @@
 
   ## Initialize the OutputFcn
   if (options.haveoutputfunction)
-    if (options.haveoutputselection)
+    if (! isempty (options.OutputSel))
       solution.retout = x(options.OutputSel,end);
     else
       solution.retout = x;
@@ -101,7 +101,7 @@
   endif
 
   ## Initialize the EventFcn
-  if (options.haveeventfunction)
+  if (! isempty (options.Events))
     ode_event_handler (options.Events, tspan(end), x,
                        "init", options.funarguments{:});
   endif
@@ -134,7 +134,7 @@
     endif
 
     err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol,
-                       options.normcontrol, x_est);
+                       strcmp (options.NormControl, "on"), x_est);
 
     ## Accept solution only if err <= 1.0
     if (err <= 1)
@@ -161,7 +161,7 @@
 
           ## Call Events function only if a valid result has been found.
           ## Stop integration if eventbreak is true.
-          if (options.haveeventfunction)
+          if (! isempty (options.Events))
             break_loop = false;
             for idenseout = 1:numel (t_caught)
               id = t_caught(idenseout);
@@ -192,7 +192,7 @@
             approxvals = interp1 ([t_old, t(t_caught), t_new],
                                   [x_old, x(:, t_caught), x_new] .',
                                   approxtime, 'linear') .';
-            if (options.haveoutputselection)
+            if (! isempty (options.OutputSel))
               approxvals = approxvals(options.OutputSel, :);
             endif
             for ii = 1:numel (approxtime)
@@ -216,7 +216,7 @@
 
         ## Call Events function only if a valid result has been found.
         ## Stop integration if eventbreak is true.
-        if (options.haveeventfunction)
+        if (! isempty (options.Events))
           solution.event = ...
           ode_event_handler (options.Events, t(istep), x(:, istep), [],
                              options.funarguments{:});
@@ -236,7 +236,7 @@
           approxvals = interp1 ([t_old, t_new],
                                 [x_old, x_new] .',
                                 approxtime, 'linear') .';
-          if (options.haveoutputselection)
+          if (! isempty (options.OutputSel))
             approxvals = approxvals(options.OutputSel, :);
           endif
           for ii = 1:numel (approxtime)