--- a/scripts/ode/module.mk	Thu Oct 06 06:54:32 2016 +0200
+++ b/scripts/ode/module.mk	Thu Oct 06 07:13:24 2016 +0200
@@ -5,12 +5,11 @@
 scripts_ode_PRIVATE_FCN_FILES = \
   scripts/ode/private/AbsRel_Norm.m \
   scripts/ode/private/integrate_adaptive.m \
-  scripts/ode/private/integrate_const.m \
-  scripts/ode/private/integrate_n_steps.m \
   scripts/ode/private/kahan.m \
   scripts/ode/private/known_option_names.m \
+  scripts/ode/private/odedefaults.m \
+  scripts/ode/private/odemergeopts.m \
   scripts/ode/private/ode_event_handler.m \
-  scripts/ode/private/ode_struct_value_check.m \
   scripts/ode/private/runge_kutta_23.m \
   scripts/ode/private/runge_kutta_45_dorpri.m \
   scripts/ode/private/runge_kutta_interpolate.m \

--- a/scripts/ode/ode23.m	Thu Oct 06 06:54:32 2016 +0200
+++ b/scripts/ode/ode23.m	Thu Oct 06 07:13:24 2016 +0200
@@ -22,7 +22,6 @@
 ## -*- texinfo -*-
 ## @deftypefn  {} {[@var{t}, @var{y}] =} ode23 (@var{fun}, @var{trange}, @var{init})
 ## @deftypefnx {} {[@var{t}, @var{y}] =} ode23 (@var{fun}, @var{trange}, @var{init}, @var{ode_opt})
-## @deftypefnx {} {[@var{t}, @var{y}] =} ode23 (@dots{}, @var{par1}, @var{par2}, @dots{})
 ## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode23 (@dots{})
 ## @deftypefnx {} {@var{solution} =} ode23 (@dots{})
 ##
@@ -40,13 +39,12 @@
 ## evaluated.  Typically, it is a two-element vector specifying the initial and
 ## final times (@code{[tinit, tfinal]}).  If there are more than two elements
 ## then the solution will also be evaluated at these intermediate time
-## instances unless the integrate function specified is
-## @code{integrate_n_steps}.
+## instances.
 ##
 ## 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}. 
+## computation may be changed by using the options @qcode{"RelTol"},
+## and @qcode{"AbsTol"},. 
 ##
 ## @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

--- a/scripts/ode/ode45.m	Thu Oct 06 06:54:32 2016 +0200
+++ b/scripts/ode/ode45.m	Thu Oct 06 07:13:24 2016 +0200
@@ -22,7 +22,6 @@
 ## -*- texinfo -*-
 ## @deftypefn  {} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init})
 ## @deftypefnx {} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init}, @var{ode_opt})
-## @deftypefnx {} {[@var{t}, @var{y}] =} ode45 (@dots{}, @var{par1}, @var{par2}, @dots{})
 ## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode45 (@dots{})
 ## @deftypefnx {} {@var{solution} =} ode45 (@dots{})
 ##
@@ -38,13 +37,12 @@
 ## evaluated.  Typically, it is a two-element vector specifying the initial and
 ## final times (@code{[tinit, tfinal]}).  If there are more than two elements
 ## then the solution will also be evaluated at these intermediate time
-## instances unless the integrate function specified is
-## @code{integrate_n_steps}.
+## instances.
 ##
 ## 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}.
+## computation may be changed by using the options @qcode{"RelTol"},
+## and @qcode{"AbsTol"}.
 ##
 ## @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

--- a/scripts/ode/private/integrate_adaptive.m	Thu Oct 06 06:54:32 2016 +0200
+++ b/scripts/ode/private/integrate_adaptive.m	Thu Oct 06 07:13:24 2016 +0200
@@ -62,7 +62,7 @@
 ##
 ## @end deftypefn
 ##
-## @seealso{integrate_const, integrate_n_steps}
+## @seealso{ode45, ode23}
 
 function solution = integrate_adaptive (stepper, order, func, tspan, x0,
                                         options)

--- a/scripts/ode/private/integrate_const.m	Thu Oct 06 06:54:32 2016 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,284 +0,0 @@
-## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {} {@var{solution} =} integrate_const (@var{@@stepper}, @var{@@func}, @var{tspan}, @var{x0}, @var{dt}, @var{options})
-##
-## This function file can be called by an ODE solver function in order to
-## integrate the set of ODEs on the interval @var{[t0,t1]} with a constant
-## timestep @var{dt}.
-##
-## The function returns a structure @var{solution} with two fieldss: @var{t}
-## and @var{y}.  @var{t} is a column vector and contains the time stamps.
-## @var{y} is a matrix in which each column refers to a different unknown
-## of the problem and the row number is the same as the @var{t} row number.
-## Thus, each row of the matrix @var{y} contains the values of all unknowns at
-## the time value contained in the corresponding row in @var{t}.
-##
-## The first input argument must be a function handle or inline function
-## representing the stepper, i.e., the function responsible for step-by-step
-## integration.  This function discriminates one method from the others.
-##
-## The second input argument is the order of the stepper.  It is needed to
-## compute the adaptive timesteps.
-##
-## The third input argument is a function handle or inline function that
-## defines the ODE:
-##
-## @ifhtml
-##
-## @example
-## @math{y' = f(t,y)}
-## @end example
-##
-## @end ifhtml
-## @ifnothtml
-## @math{y' = f(t,y)}.
-## @end ifnothtml
-##
-## The fourth input argument is the time vector which defines the integration
-## interval, i.e., @var{[tspan(1), tspan(end)]} and all intermediate elements
-## are taken as times at which the solution is required.
-##
-## The fourth argument contains the initial conditions for the ODEs.
-##
-## The fifth input argument represents the fixed timestep and the last input
-## argument contains some options that may be needed for the stepper.
-## @end deftypefn
-##
-## @seealso{integrate_adaptive, integrate_n_steps}
-
-function solution = integrate_const (stepper, func, tspan, x0, dt, options)
-
-  solution = struct ();
-
-  ## first values for time and solution
-  t = tspan(1);
-  x = x0(:);
-
-  direction = odeget (options, "direction", [], "fast");
-  if (sign (dt) != direction)
-    error ("Octave:invalid-input-arg",
-           "option 'InitialStep' has a wrong sign");
-  endif
-
-  ## setting parameters
-  k = length (tspan);
-  counter = 2;
-  comp = 0.0;
-  tk = tspan(1);
-  options.comp = comp;
-
-  ## Initialize the OutputFcn
-  if (options.haveoutputfunction)
-    if (options.haveoutputselection)
-      solution.retout = x(options.OutputSel,end);
-    else
-      solution.retout = x;
-    endif
-    feval (options.OutputFcn, tspan, solution.retout, "init",
-           options.funarguments{:});
-  endif
-
-  ## Initialize the EventFcn
-  if (options.haveeventfunction)
-    ode_event_handler (options.Events, t(end), x, "init",
-                         options.funarguments{:});
-  endif
-
-  solution.cntloop = 2;
-  solution.cntcycles = 1;
-  cntiter = 0;
-  solution.unhandledtermination = true;
-  solution.cntsave = 2;
-
-  z = t;
-  u = x;
-  k_vals = feval (func, t , x, options.funarguments{:});
-
-  while (counter <= k)
-    ## computing the integration step from t to t+dt
-    [s, y, ~, k_vals] = stepper (func, z(end), u(:,end), dt, options, k_vals);
-
-    [tk, comp] = kahan (tk,comp, dt);
-    options.comp = comp;
-    s(end) = tk;
-
-    if (options.havenonnegative)
-      x(options.NonNegative,end) = abs (x(options.NonNegative,end));
-      y(options.NonNegative,end) = abs (y(options.NonNegative,end));
-      y_est(options.NonNegative,end) = abs (y_est(options.NonNegative,end));
-    endif
-
-    if (options.haveoutputfunction && options.haverefine)
-      SaveVUForRefine = u(:,end);
-    endif
-
-    ## values on this interval for time and solution
-    z = [t(end);s];
-    u = [x(:,end),y];
-
-    ## if next tspan value is caught, update counter
-    if ((z(end) == tspan(counter))
-        || (abs (z(end) - tspan(counter)) /
-            (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) )
-      counter += 1;
-
-    ## if there is an element in time vector at which the solution is required
-    ## the program must compute this solution before going on with next steps
-    elseif (direction * z(end) > direction * tspan(counter) )
-      ## initializing counter for the following cycle
-      i = 2;
-      while (i <= length (z))
-
-        ## if next tspan value is caught, update counter
-        if ((counter <= k)
-            && (((z(i) == tspan(counter))
-                 || (abs (z(i) - tspan(counter)) /
-                     (max (abs (z(i)), abs (tspan(counter)))) < 8*eps))) )
-          counter += 1;
-        endif
-        ## else, loop until there are requested values inside this subinterval
-        while ((counter <= k)
-               && direction * z(i) > direction * tspan(counter) )
-          ## add the interpolated value of the solution
-          u = [u(:,1:i-1),u(:,i-1) + (tspan(counter)-z(i-1))/(z(i)-z(i-1))* ...
-              (u(:,i)-u(:,i-1)),u(:,i:end)];
-          ## add the time requested
-          z = [z(1:i-1);tspan(counter);z(i:end)];
-
-          ## update counters
-          counter += 1;
-          i += 1;
-        endwhile
-
-        ## if new time requested is not out of this interval
-        if (counter <= k && direction * z(end) > direction * tspan(counter))
-          ## update the counter
-          i += 1;
-        else
-          ## else, stop the cycle and go on with the next iteration
-          i = length (z)+1;
-        endif
-
-      endwhile
-    endif
-
-    x = [x,u(:,2:end)];
-    t = [t;z(2:end)];
-    solution.cntsave += 1;
-    solution.cntloop += 1;
-    cntiter = 0;
-
-    ## Call OutputFcn only if a valid result has been found.
-    ## Stop integration if function returns false.
-    if (options.haveoutputfunction)
-      for cnt = 0:options.Refine  # Approximation between told and t
-        if (options.haverefine)   # Do interpolation
-          approxtime = (cnt + 1) / (options.Refine + 2);
-          approxvals = (1 - approxtime) * SaveVUForRefine ...
-                        + (approxtime) * y(:,end);
-          approxtime = s(end) + approxtime*dt;
-        else
-          approxvals = x(:,end);
-          approxtime = t(end);
-        endif
-        if (options.haveoutputselection)
-          approxvals = approxvals(options.OutputSel);
-        endif
-        pltret = feval (options.OutputFcn, approxtime, approxvals, [],
-                         options.funarguments{:});
-        if (pltret)  # Leave refinement loop
-          break;
-        endif
-      endfor
-      if (pltret)  # Leave main loop
-        solution.unhandledtermination = false;
-        break;
-      endif
-    endif
-
-    ## Call Events function only if a valid result has been found.
-    ## Stop integration if eventbreak is true.
-    if (options.haveeventfunction)
-      solution.event = ode_event_handler (options.Events, t(end), x(:,end),
-                                             [], options.funarguments{:});
-      if (! isempty (solution.event{1}) && solution.event{1} == 1)
-        t(solution.cntloop-1,:) = solution.event{3}(end,:);
-        x(:,solution.cntloop-1) = solution.event{4}(end,:)';
-        solution.unhandledtermination = false;
-        break;
-      endif
-    endif
-
-    ## Update counters that count the number of iteration cycles
-    solution.cntcycles += 1;  # Needed for cost statistics
-    cntiter += 1;             # Needed to find iteration problems
-
-    ## Stop solving because, in the last 5,000 steps, no successful valid
-    ## value has been found
-    if (cntiter >= 5_000)
-      error (["integrate_const: Solving was not successful. ", ...
-              " The iterative integration loop exited at time", ...
-              " t = %f before the endpoint at tend = %f was reached. ", ...
-              " This happened because the iterative integration loop", ...
-              " did not find a valid solution at this time stamp. ", ...
-              " Try to reduce the value of 'InitialStep' and/or", ...
-              " 'MaxStep' with the command 'odeset'."],
-             s(end), tspan(end));
-    endif
-
-    ## if this is the last iteration, save the length of last interval
-    if (counter > k)
-      j = length (z);
-    endif
-
-  endwhile
-
-  ## Check if integration of the ode has been successful
-  if (direction * z(end) < direction * tspan(end))
-    if (solution.unhandledtermination == true)
-      error ("integrate_const:unexpected_termination",
-             [" Solving was not successful. ", ...
-              " The iterative integration loop exited at time", ...
-              " t = %f before the endpoint at tend = %f was reached. ", ...
-              " This may happen if the stepsize becomes too small. ", ...
-              " Try to reduce the value of 'InitialStep'", ...
-              " and/or 'MaxStep' with the command 'odeset'."],
-             z(end), tspan(end));
-    else
-      warning ("integrate_const:unexpected_termination",
-               ["Solver was stopped by a call of 'break'", ...
-                " in the main iteration loop at time", ...
-                " t = %f before the endpoint at tend = %f was reached. ", ...
-                " This may happen because the @odeplot function", ...
-                " returned 'true' or the @event function returned 'true'."],
-               z(end), tspan(end));
-    endif
-  endif
-
-  ## compute how many values are out of time inerval
-  d = direction * t((end-(j-1)):end) > direction * tspan(end) * ones (j, 1);
-  f = sum (d);
-
-  ## remove not-requested values of time and solution
-  solution.t = t(1:end-f);
-  solution.x = x(:,1:end-f)';
-
-endfunction
-

--- a/scripts/ode/private/integrate_n_steps.m	Thu Oct 06 06:54:32 2016 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,223 +0,0 @@
-## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {} {@var{solution} =} integrate_n_steps (@var{@@stepper}, @var{@@func}, @var{t0}, @var{x0}, @var{dt}, @var{n}, @var{options})
-##
-## This function file can be called by an ODE solver function in order to
-## integrate the set of ODEs on the interval @var{[t0,t0 + n*dt]} with a
-## constant timestep dt and on a fixed number of steps.
-##
-## The function returns a structure @var{solution} with two fieldss: @var{t}
-## and @var{y}.  @var{t} is a column vector and contains the time stamps.
-## @var{y} is a matrix in which each column refers to a different unknown
-## of the problem and the row number is the same as the @var{t} row number.
-## Thus, each row of the matrix @var{y} contains the values of all unknowns at
-## the time value contained in the corresponding row in @var{t}.
-##
-## The first input argument must be a function handle or inline function
-## representing the stepper, i.e., the function responsible for step-by-step
-## integration.  This function discriminates one method from the others.
-##
-## The second input argument is the order of the stepper.  It is needed to
-## compute the adaptive timesteps.
-##
-## The third input argument is a function handle or inline function that
-## defines the ODE:
-##
-## @ifhtml
-##
-## @example
-## @math{y' = f(t,y)}
-## @end example
-##
-## @end ifhtml
-## @ifnothtml
-## @math{y' = f(t,y)}.
-## @end ifnothtml
-##
-## The third input argument is the starting point for the integration.
-##
-## The fourth argument contains the initial conditions for the ODEs.
-##
-## The fifth input argument represents the fixed timestep while the sixth
-## contains the number of integration steps.
-##
-## The last argument is a struct with the options that may be needed by the
-## stepper.
-## @end deftypefn
-##
-## @seealso{integrate_adaptive, integrate_const}
-
-function solution = integrate_n_steps (stepper, func, t0, x0, dt, n, options)
-
-  solution = struct ();
-
-  ## first values for time and solution
-  x = x0(:);
-  t = t0;
-
-  direction = odeget (options, "direction", [], "fast");
-  if (sign (dt) != direction)
-    error ("Octave:invalid-input-arg", "option 'InitialStep' has a wrong sign");
-  endif
-
-  comp = 0.0;
-  tk = t0;
-  options.comp = comp;
-
-  ## Initialize the OutputFcn
-  if (options.haveoutputfunction)
-    if (options.haveoutputselection)
-      solution.retout = x(options.OutputSel,end);
-    else
-      solution.retout = x;
-    endif
-    feval (options.OutputFcn, tspan, solution.retout, "init",
-           options.funarguments{:});
-  endif
-
-  ## Initialize the EventFcn
-  if (options.haveeventfunction)
-    ode_event_handler (options.Events, t(end), x, "init",
-                         options.funarguments{:});
-  endif
-
-  solution.cntloop = 2;
-  solution.cntcycles = 1;
-  cntiter = 0;
-  solution.unhandledtermination = true;
-  solution.cntsave = 2;
-
-  z = t;
-  u = x;
-  k_vals = feval (func, t , x, options.funarguments{:});
-
-  for i = 1:n
-    ## Compute the integration step from t to t+dt
-    [s, y, ~, k_vals] = stepper (func, z(end), u(:,end), dt, options, k_vals);
-
-    [tk, comp] = kahan (tk, comp, dt);
-    options.comp = comp;
-    s(end) = tk;
-
-    if (options.havenonnegative)
-      x(options.NonNegative,end) = abs (x(options.NonNegative,end));
-      y(options.NonNegative,end) = abs (y(options.NonNegative,end));
-    endif
-
-    if (options.haveoutputfunction && options.haverefine)
-      SaveVUForRefine = u(:,end);
-    endif
-
-    ## values on this interval for time and solution
-    z = [t(end);s];
-    u = [x(:,end),y];
-
-    x = [x,u(:,2:end)];
-    t = [t;z(2:end)];
-    solution.cntsave += 1;
-    solution.cntloop += 1;
-    cntiter = 0;
-
-    ## Call OutputFcn only if a valid result has been found.
-    ## Stop integration if function returns false.
-    if (options.haveoutputfunction)
-      for cnt = 0:options.Refine  # Approximation between told and t
-        if (options.haverefine)   # Do interpolation
-          approxtime = (cnt + 1) / (options.Refine + 2);
-          approxvals = (1 - approxtime) * SaveVUForRefine ...
-                        + (approxtime) * y(:,end);
-          approxtime = s(end) + approxtime*dt;
-        else
-          approxvals = x(:,end);
-          approxtime = t(end);
-        endif
-        if (options.haveoutputselection)
-          approxvals = approxvals(options.OutputSel);
-        endif
-        pltret = feval (options.OutputFcn, approxtime, approxvals, [],
-                        options.funarguments{:});
-        if (pltret)  # Leave refinement loop
-          break;
-        endif
-      endfor
-      if (pltret)  # Leave main loop
-        solution.unhandledtermination = false;
-        break;
-      endif
-    endif
-
-    ## Call Events function only if a valid result has been found.
-    ## Stop integration if eventbreak is true.
-    if (options.haveeventfunction)
-      solution.event = ode_event_handler (options.Events, t(end), x(:,end),
-                                             [], options.funarguments{:});
-      if (! isempty (solution.event{1}) && solution.event{1} == 1)
-        t(solution.cntloop-1,:) = solution.event{3}(end,:);
-        x(:,solution.cntloop-1) = solution.event{4}(end,:)';
-        solution.unhandledtermination = false;
-        break;
-      endif
-    endif
-
-    ## Update counters that count the number of iteration cycles
-    solution.cntcycles += 1;  # Needed for cost statistics
-    cntiter += 1;             # Needed to find iteration problems
-
-    ## Stop solving because, in the last 5,000 steps, no successful valid
-    ## value has been found
-    if (cntiter >= 5_000)
-      error (["integrate_n_steps: Solving was not successful. ", ...
-              " The iterative integration loop exited at time", ...
-              " t = %f before the endpoint at tend = %f was reached. ", ...
-              " This happened because the iterative integration loop", ...
-              " did not find a valid solution at this time stamp. ", ...
-              " Try to reduce the value of 'InitialStep' and/or", ...
-              " 'MaxStep' with the command 'odeset'."],
-             s(end), tspan(end));
-    endif
-  endfor
-
-  ## Check if integration of the ode has been successful
-  #if (direction * z(end) < direction * tspan(end))
-  #  if (solution.unhandledtermination == true)
-  #   error ("integrate_n_steps:unexpected_termination",
-  #          [" Solving was not successful. ", ...
-  #           " The iterative integration loop exited at time", ...
-  #           " t = %f before the endpoint at tend = %f was reached. ", ...
-  #           " This may happen if the stepsize becomes too small. ", ...
-  #           " Try to reduce the value of 'InitialStep'", ...
-  #           " and/or 'MaxStep' with the command 'odeset'."],
-  #           z(end), tspan(end));
-  #  else
-  #   warning ("integrate_n_steps:unexpected_termination",
-  #            ["Solver was stopped by a call of 'break'", ...
-  #             " in the main iteration loop at time", ...
-  #             " t = %f before the endpoint at tend = %f was reached. ", ...
-  #             " This may happen because the @odeplot function", ...
-  #             " returned 'true' or the @event function returned 'true'."],
-  #             z(end), tspan(end));
-  #  endif
-  #endif
-
-  solution.t = t;
-  solution.x = x';
-
-endfunction
-

--- a/scripts/ode/private/ode_struct_value_check.m	Thu Oct 06 06:54:32 2016 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,294 +0,0 @@
-## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
-## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net>
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn  {} {} ode_struct_value_check (@var{"caller"}, @var{ode_struct})
-## @deftypefnx {} {} ode_struct_value_check (@var{"caller"), @var{ode_struct}, @var{"solver"})
-## @deftypefnx {} {@var{ode_struct} =} ode_struct_value_check (@dots{})
-##
-## Validate the fields and values in the ODE options structure
-## @var{ode_struct}.
-##
-## The first argument @var{caller} is a string with the name of the calling
-## function so that warning and error messages properly display the source
-## of any problems.
-##
-## The second argument @var{ode_struct} is a structure with fields and values
-## that configure the ODE solvers (@pxref{XREFodeset,,odeset}).
-##
-## The optional third argument @var{"solver"} is a string with the name of a
-## specific ODE solver.  This extra information can enable more extensive value
-## validation for certain options.
-##
-## The function does not modify any of the field names or field values, but
-## terminates with an error if an invalid value is found.
-##
-## Normally the function is called with no output.  However, the input struct
-## is passed unmodified to the output for certain solvers which expect to
-## receive the validated ODE structure returned.
-## @end deftypefn
-##
-## @seealso{odeset, odeget}
-
-function ode_struct = ode_struct_value_check (caller, ode_struct, solver = "")
-
-  for [val, opt] = ode_struct  # Cycle over all fields
-
-    switch (opt)
-
-      case "AbsTol"
-        if (! isempty (val))
-          if (! isnumeric (val) || ! isreal (val)
-              || ! isvector (val) || any (val <= 0))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "BDF"
-        if (! isempty (val))
-          if (! strcmp (val, "on") && ! strcmp (val, "off"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "Events"
-        if (! isempty (val))
-          if (! isa (val, "function_handle"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "InitialSlope"
-        if (! isempty (val))
-          if (! ischar (val)
-              && (! isnumeric (val) || (! isvector (val) && ! isreal (val))))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "InitialStep"
-        if (! isempty (val))
-          if (! isnumeric (val) || ! isreal (val) || ! isscalar (val)
-              || val <= 0)
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "Jacobian"
-        if (! isempty (val))
-          if (! isnumeric (val))
-            if (! isa (val, "function_handle") && ! iscell (val))
-              error ("Octave:invalid-input-arg",
-                     [caller ": invalid value assigned to field '%s'"], opt);
-            endif
-          endif
-        endif
-
-      case "JConstant"
-        if (! isempty (val))
-          if (! strcmp (val, "on") && ! strcmp (val, "off"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "JPattern"
-        if (! isempty (val))
-          if (! isnumeric (val) && ! isvector (val))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "Mass"
-        if (! isempty (val))
-          if ((! isnumeric (val) || ! ismatrix (val))
-              && ! isa (val, "function_handle"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "MassConstant"
-        if (! isempty (val))
-          if (! strcmp (val, "on") && ! strcmp (val, "off"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "MassSingular"
-        if (! isempty (val))
-          if (! any (strcmp (val, {"yes", "no", "maybe"})))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "MaxOrder"
-        if (! isempty (val))
-          if (! isnumeric (val)
-              || val != fix (val) || val <= 0 || val >= 8)
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "MaxStep"
-        if (! isempty (val))
-          if (! isnumeric (val) || ! isscalar (val) || val <= 0)
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "MStateDependence"
-        if (! isempty (val))
-          if (! any (strcmp (val, {"none", "weak", "strong"})))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "MvPattern"
-        if (! isempty (val))
-          if (! isnumeric (val) && ! isvector (val))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "NonNegative"
-        if (! isempty (val))
-          if (! isnumeric (val) || ! isvector (val)
-              || any (val <= 0) || any (val != fix (val)))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "NormControl"
-        if (! isempty (val))
-          if (! strcmp (val, "on") && ! strcmp (val, "off"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "OutputFcn"
-        if (! isempty (val))
-          if (! isa (val, "function_handle"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "OutputSel"
-        if (! isempty (val))
-          if (! isnumeric (val) || ! isvector (val))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "Refine"
-        if (! isempty (val))
-          if (! isnumeric (val) || ! isscalar (val)
-              || val != fix (val)  || val < 0 || val > 5)
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "RelTol"
-        if (! isempty (val))
-          if (! isnumeric (val) || ! isreal (val) || any (val <= 0))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-          if (any (strcmp (solver, {"ode23", "ode23d", "ode45", "ode45d",
-                                    "ode54", "ode54d", "ode78", "ode78d"})))
-            if (! isscalar (val))
-              error ("Octave:invalid-input-arg",
-                     [caller ": invalid value assigned to field '%s'"], opt);
-            endif
-          endif
-        endif
-
-      case "Stats"
-        if (! isempty (val))
-          if (! strcmp (val, "on") && ! strcmp (val, "off"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "Vectorized"
-        if (! isempty (val))
-          if (! strcmp (val, "on") && ! strcmp (val, "off"))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      case "TimeStepSize"
-        if (! isempty (val))
-         if (! isscalar (val))
-              error ("Octave:invalid-input-arg",
-                     [caller ": invalid value assigned to field '%s'"], opt);
-         endif
-        endif
-
-      case "TimeStepNumber"
-        if (! isempty (val))
-          if (! isscalar (val))
-            error ("Octave:invalid-input-arg",
-                   [caller ": invalid value assigned to field '%s'"], opt);
-          endif
-        endif
-
-      otherwise
-        warning ("Octave:invalid-input-arg",
-                 [caller ": unknown field '%s' in ODE options\n"], opt);
-    endswitch
-  endfor
-
-endfunction
-
-
-%!demo
-%! # Return the checked ODE options structure that is created by
-%! # the command odeset.
-%!
-%! ode_struct_value_check (odeset);
-
-%!demo
-%! # Create the ODE options structure A with odeset and check it
-%! # with ode_struct_value_check.  This actually is unnecessary
-%! # because odeset automatically calls ode_struct_value_check before
-%! # returning.
-%!
-%! A = odeset ();
-%! ode_struct_value_check (A);
-

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/ode/private/odedefaults.m	Thu Oct 06 07:13:24 2016 +0200
@@ -0,0 +1,96 @@
+## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+function [defaults, classes, attributes] = odedefaults (n, t0, tf)
+
+defaults = odeset ('AbsTol', 1e-6,
+                   'BDF', 'off',
+                   'Events', [],
+                   'InitialSlope', zeros(n,1),
+                   'InitialStep', [],
+                   'Jacobian', [],
+                   'JConstant', 'off',
+                   'JPattern', [],
+                   'Mass', [],
+                   'MassConstant', 'off',
+                   'MassSingular', 'maybe',
+                   'MaxOrder', 5,
+                   'MaxStep', 0.1*abs(t0-tf),
+                   'MStateDependence', 'weak',
+                   'MvPattern', [],
+                   'NonNegative', [],
+                   'NormControl', 'off',
+                   'OutputFcn', [],
+                   'OutputSel', [],
+                   'Refine', 1,  
+                   'RelTol', 1e-3,
+                   'Stats', 'off',
+                   'Vectorized', 'off');
+ 
+classes = odeset ('Abstol', {"float"},
+                  'BDF', "char",
+                  'Events', {"function_handle"},
+                  'InitialSlope', {"float"},
+                  'InitialStep', {"float"},
+                  'Jacobian', {"float", "function_handle", "cell"},
+                  'JConstant', "char",
+                  'JPattern', {"float"},
+                  'Mass', {"float", "function_handle"},
+                  'MassConstant', "char",
+                  'MassSingular', "char",
+                  'MaxOrder', {"float"},
+                  'MaxStep', {"float"},
+                  'MStateDependence', "char",
+                  'MvPattern', {"float"},
+                  'NonNegative', {"float"},
+                  'NormControl', "char",
+                  'OutputFcn', {"function_handle"},
+                  'OutputSel', {"float"},
+                  'Refine', {"float"},
+                  'RelTol', {"float"},
+                  'Stats', "char",
+                  'Vectorized', "char");
+
+
+##FIXME: How can I check Jacobian where it's a cell????? Maybe it's better to check it inside the solver
+##FIXME: Vectorized can be a cell of stings
+attributes = odeset ('AbsTol', {"real", "vector", "positive"},
+                     'BDF', {"on", "off"},
+                     'Events', {},
+                     'InitialSlope', {"real", "vector", "numel", n},
+                     'InitialStep', {"positive", "scalar"},
+                     'Jacobian', {},
+                     'JConstant', {"on", "off"},
+                     'JPattern', {"vector"},
+                     'Mass', {},
+                     'MassConstant', {"on", "off"},
+                     'MassSingular', {"no", "maybe", "yes"},
+                     'MaxOrder', {">=", 0, "<=", 5, "integer"},
+                     'MaxStep', {"positive", "scalar", "real"},
+                     'MStateDependence', {"weak", "strong", "none"},
+                     'MvPattern', {"vector"},
+                     'NonNegative', {"vector", "integer", "positive"},
+                     'NormControl', {"on", "off"},
+                     'OutputFcn', {},
+                     'OutputSel', {"vector", "integer", "positive",...
+                                   ">", 0, "<=", n},
+                     'Refine', {"scalar", ">", 0, "integer"},
+                     'RelTol', {"scalar", "positive", "real"},
+                     'Stats', {"on", "off"},
+                     'Vectorized', {"on", "off"});
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/ode/private/odemergeopts.m	Thu Oct 06 07:13:24 2016 +0200
@@ -0,0 +1,42 @@
+## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+function options = odemergeopts  (useroptions, options, classes,
+                                  attributes, fun_name);
+
+  for [value, key] = options;
+
+    if (isfield (useroptions, key) && ! isempty (useroptions.(key)))
+
+      if (! strcmp (classes.(key), "char"))
+        validateattributes (useroptions.(key), classes.(key),
+                            attributes.(key), fun_name, key);
+
+      elseif (ischar (useroptions.(key)))
+        validatestring (useroptions.(key), attributes.(key), fun_name, key);
+
+      else
+        error ("Octave:invalid-input-arg",
+                [fun_name ": invalid value assigned to field '%s'"], key);
+      endif
+      
+    options.(key) = useroptions.(key);
+    
+    endif
+  endfor
+endfunction