Mercurial > octave

changeset 22649:2d9393be9a18

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Periodic merge of stable to default.
author Rik <rik@octave.org>
date Wed, 19 Oct 2016 15:20:33 -0700
parents 859d48b5648d (current diff) 9bc03a3f7a34 (diff)
children 65d84345abfc
files
diffstat 8 files changed, 210 insertions(+), 221 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ode/ode23.m	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/ode23.m	Wed Oct 19 15:20:33 2016 -0700
@@ -89,15 +89,6 @@
 ## @seealso{odeset, odeget, ode45}
 ## @end deftypefn
 
-## FIXME: We store ChangeLog information in Mercurial.  Can this be deleted?
-## ChangeLog:
-##   20010703 the function file "ode23.m" was written by Marc Compere
-##     under the GPL for the use with this software.  This function has been
-##     taken as a base for the following implementation.
-##   20060810, Thomas Treichl
-##     This function was adapted to the new syntax that is used by the
-##     new OdePkg for Octave and is compatible to Matlab's ode23.
-
 function varargout = ode23 (fun, trange, init, varargin)
 
   if (nargin < 3)
@@ -161,12 +152,6 @@
 
   ## Start preprocessing, have a look which options are set in odeopts,
   ## check if an invalid or unused option is set.
-  ## FIXME: Why persistent?  Won't these have different values for every
-  ##        run of ode23?
-  persistent defaults   = [];
-  persistent classes    = [];
-  persistent attributes = [];
-
   [defaults, classes, attributes] = odedefaults (numel (init),
                                                  trange(1), trange(end));
 
@@ -243,8 +228,7 @@
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, solution.t(end),
-           solution.x(end,:)', "done", odeopts.funarguments{:});
+    feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
     ode_event_handler (odeopts.Events, solution.t(end),
@@ -312,9 +296,6 @@
 %!   err(i) = norm (y .* exp (t) - 1, Inf);
 %! endfor
 %!
-%! ## Estimate order numerically
-%! p = diff (log (err)) ./ diff (log (h))
-%!
 %! ## Estimate order visually
 %! loglog (h, tol, "-ob",
 %!         h, err, "-b",
@@ -326,6 +307,9 @@
 %! title ("Convergence plot for ode23");
 %! legend ("imposed tolerance", "ode23 (relative) error",
 %!         "order 2", "order 3", "location", "northwest");
+%!
+%! ## Estimate order numerically
+%! p = diff (log (err)) ./ diff (log (h))
 
 ## We are using the "Van der Pol" implementation for all tests that are done
 ## for this function.
@@ -360,15 +344,21 @@
 %!  mas = sparse ([1, 0; 0, 1]);    # A sparse dummy matrix
 %!endfunction
 %!function out = fout (t, y, flag, varargin)
-%!  if (regexp (char (flag), "init") == 1)
-%!    if (any (size (t) != [2, 1])) error ('"fout" step "init"'); endif
+%!  out = false;
+%!  if (strcmp (flag, "init"))
+%!    if (! isequal (size (t), [2, 1]))
+%!      error ('fout: step "init"');
+%!    endif
 %!  elseif (isempty (flag))
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "calc"'); endif
-%!    out = false;
-%!  elseif (regexp (char (flag), "done") == 1)
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "done"'); endif
+%!    if (! isequal (size (t), [1, 1]))
+%!      error ('fout: step "calc"');
+%!    endif
+%!  elseif (strcmp (flag, "done"))
+%!    if (! isempty (t))
+%!      warning ('fout: step "done"');
+%!    endif
 %!  else
-%!    error ('"fout" invalid flag');
+%!    error ("fout: invalid flag <%s>", flag);
 %!  endif
 %!endfunction
 %!
@@ -376,8 +366,8 @@
 %! [t, y] = ode23 (@fpol, [0 2], [2 0]);
 %! assert ([t(end), y(end,:)], [2, fref], 1e-3);
 %!test  # anonymous function instead of real function
-%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
-%! [t, y] = ode23 (fvdb, [0 2], [2 0]);
+%! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
+%! [t, y] = ode23 (fvdp, [0 2], [2 0]);
 %! assert ([t(end), y(end,:)], [2, fref], 1e-3);
 %!test  # extra input arguments passed through
 %! [t, y] = ode23 (@fpol, [0 2], [2 0], 12, 13, "KL");
@@ -484,9 +474,9 @@
 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
 
 ## FIXME: Missing tests.
-## 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 input validation
 %!error ode23 ()
--- a/scripts/ode/ode45.m	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/ode45.m	Wed Oct 19 15:20:33 2016 -0700
@@ -143,11 +143,6 @@
 
   ## Start preprocessing, have a look which options are set in odeopts,
   ## check if an invalid or unused option is set
-  ## FIXME: Why persistent when it is changed with every run of ode45?
-  persistent defaults   = [];
-  persistent classes    = [];
-  persistent attributes = [];
-
   [defaults, classes, attributes] = odedefaults (numel (init),
                                                  trange(1), trange(end));
 
@@ -224,8 +219,7 @@
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, solution.t(end),
-           solution.x(end,:)', "done", odeopts.funarguments{:});
+    feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
     ode_event_handler (odeopts.Events, solution.t(end),
@@ -293,9 +287,6 @@
 %!   err(i) = norm (y .* exp (t) - 1, Inf);
 %! endfor
 %!
-%! ## Estimate order numerically
-%! p = diff (log (err)) ./ diff (log (h))
-%!
 %! ## Estimate order visually
 %! loglog (h, tol, "-ob",
 %!         h, err, "-b",
@@ -307,6 +298,9 @@
 %! title ("Convergence plot for ode45");
 %! legend ("imposed tolerance", "ode45 (relative) error",
 %!         "order 4", "order 5", "location", "northwest");
+%!
+%! ## Estimate order numerically
+%! p = diff (log (err)) ./ diff (log (h))
 
 ## We are using the Van der Pol equation for all tests that are done
 ## for this function.
@@ -341,15 +335,21 @@
 %!  mas = sparse ([1, 0; 0, 1]);   # A sparse dummy matrix
 %!endfunction
 %!function out = fout (t, y, flag, varargin)
-%!  if (regexp (char (flag), "init") == 1)
-%!    if (any (size (t) != [2, 1])) error ('"fout" step "init"'); endif
+%!  out = false;
+%!  if (strcmp (flag, "init"))
+%!    if (! isequal (size (t), [2, 1]))
+%!      error ('fout: step "init"');
+%!    endif
 %!  elseif (isempty (flag))
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "calc"'); endif
-%!    out = false;
-%!  elseif (regexp (char (flag), 'done') == 1)
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "done"'); endif
+%!    if (! isequal (size (t), [1, 1]))
+%!      error ('fout: step "calc"');
+%!    endif
+%!  elseif (strcmp (flag, "done"))
+%!    if (! isempty (t))
+%!      warning ('fout: step "done"');
+%!    endif
 %!  else
-%!    error ('"fout" invalid flag');
+%!    error ("fout: invalid flag <%s>", flag);
 %!  endif
 %!endfunction
 %!
@@ -360,8 +360,8 @@
 %! [t, y] = ode45 (@fpol, [0 2], [2 0]);
 %! assert (size (t) < 20);
 %!test  # anonymous function instead of real function
-%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
-%! [t, y] = ode45 (fvdb, [0 2], [2 0]);
+%! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
+%! [t, y] = ode45 (fvdp, [0 2], [2 0]);
 %! assert ([t(end), y(end,:)], [2, fref], 1e-2);
 %!test  # string instead of function
 %! [t, y] = ode45 ("fpol", [0 2], [2 0]);
@@ -482,9 +482,9 @@
 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3);
 
 ## FIXME: Missing tests.
-## test for MvPattern option is missing
 ## test for InitialSlope option is missing
 ## test for MaxOrder option is missing
+## test for MvPattern option is missing
 
 %!error ode45 ()
 %!error ode45 (1)
--- a/scripts/ode/odeget.m	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/odeget.m	Wed Oct 19 15:20:33 2016 -0700
@@ -37,9 +37,7 @@
 ## @seealso{odeset}
 ## @end deftypefn
 
-## FIXME: 4th input argument "opt" is undocumented.
-
-function val = odeget (ode_opt, field, default = [], opt = "")
+function val = odeget (ode_opt, field, default = [])
 
   validateattributes (ode_opt, {"struct"}, {"nonempty"});
   validateattributes (field, {"char"}, {"nonempty"});
@@ -58,9 +56,9 @@
 
 
 %!demo
-%! # Return the manually changed value RelTol of the ODE options
-%! # structure A.  If RelTol wouldn't have been changed then an
-%! # empty matrix value would have been returned.
+%! ## Return the manually changed value RelTol of the ODE options
+%! ## structure 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", [])
@@ -73,7 +71,9 @@
 %!assert (odeget (odeset (), "Mass"), [])
 %!assert (odeget (odeset (), "AbsTol", 1e-9), 1e-9)
 %!assert (odeget (odeset ("AbsTol", 1e-9), "AbsTol", []), 1e-9)
-%!assert (odeget (odeset ("foo", 42), "foo"), 42)
+%!test
+%! warning ("off", "Octave:invalid-input-arg", "local");
+%! assert (odeget (odeset ("foo", 42), "foo"), 42);
 
 %!error odeget ()
 %!error odeget (1)
--- a/scripts/ode/odeplot.m	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/odeplot.m	Wed Oct 19 15:20:33 2016 -0700
@@ -19,80 +19,93 @@
 ## Author: Thomas Treichl <treichl@users.sourceforge.net>
 
 ## -*- texinfo -*-
-## @deftypefn {} {@var{ret} =} odeplot (@var{t}, @var{y}, @var{flag})
+## @deftypefn {} {@var{stop_solve} =} odeplot (@var{t}, @var{y}, @var{flag})
 ##
-## Open a new figure window and plot input @var{y} over time during the
-## solving of an ode problem.
+## Open a new figure window and plot the solution of an ode problem at each
+## time step during the integration.
 ##
-## The input @var{y} is a column vector.  The types and values of the input
-## parameter @var{t} and output parameter @var{ret} depend on the input
-## @var{flag} that is of type string.  If @var{flag} is
+## The types and values of the input parameters @var{t} and @var{y} depend on
+## the input @var{flag} that is of type string.  Valid values of @var{flag}
+## are:
 ##
 ## @table @option
 ## @item @qcode{"init"}
-## then @var{t} must be a column vector of length 2 with the first and
-## the last time step;  Nothing is returned from this function.
+## The input @var{t} must be a column vector of length 2 with the first and
+## last time step (@code{[@var{tfirst} @var{tlast}]}.  The input @var{y}
+## contains the initial conditions for the ode problem (@var{y0}).
 ##
 ## @item @qcode{""}
-## then @var{t} must be a scalar double specifying the actual time step;
-## The return value is false (resp. value 0) for @qcode{"not stop solving"}.
+## The input @var{t} must be a scalar double specifying the time for which
+## the solution in input @var{y} was calculated.
 ##
 ## @item @qcode{"done"}
-## then @var{t} must be a scalar double specifying the last time step;
-## Nothing is returned from this function.
+## The inputs should be empty, but are ignored if they are present.
 ## @end table
 ##
-## This function is called by an ode solver function if it was specified in
-## an options structure with @code{odeset}.  This function is an internal
-## helper function; It should never be necessary for this function to be
-## directly called by a user.  There is only minimal error detection
-## implemented in order to to achieve the highest performance.
+## @code{odeplot} always returns false, i.e., don't stop the ode solver.
 ##
-## For example, solve an anonymous implementation of the
+## Example: solve an anonymous implementation of the
 ## @nospell{@qcode{"Van der Pol"}} equation and display the results while
-## solving
+## solving.
 ##
 ## @example
 ## @group
-## fvdb = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
+## fvdp = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
 ##
 ## opt = odeset ("OutputFcn", @@odeplot, "RelTol", 1e-6);
-## sol = ode45 (fvdb, [0 20], [2 0], opt);
+## sol = ode45 (fvdp, [0 20], [2 0], opt);
 ## @end group
 ## @end example
 ##
-## @seealso{odeset, odeget}
+## Background Information:
+## This function is called by an ode solver function if it was specified in
+## the @qcode{"OutputFcn"} property of an options structure created with
+## @code{odeset}.  The ode solver will initially call the function with the
+## syntax @code{odeplot ([@var{tfirst}, @var{tlast}], @var{y0}, "init")}.  The
+## function initializes internal variables, creates a new figure window, and
+## sets the x limits of the plot.  Subsequently, at each time step during the
+## integration the ode solver calls @code{odeplot (@var{t}, @var{y}, [])}.
+## At the end of the solution the ode solver calls
+## @code{odeplot ([], [], "end")} so that odeplot can perform any clean-up
+## actions required.
+## @seealso{odeset, odeget, ode23, ode45}
 ## @end deftypefn
 
-function ret = odeplot (t, y, flag, varargin)
+function stop_solve = odeplot (t, y, flag)
 
-  ## No input argument check is done for a higher processing speed
-  persistent fig told yold counter;
+  ## No input argument checking is done for better performance
+  persistent hlines num_lines told yold;
+  persistent idx = 1;   # Don't remove.  Required for Octave parser.
+
+  ## odeplot never stops the integration
+  stop_solve = false;
 
-  if (strcmp (flag, "init"))
-    ## Nothing to return, t is either the time slot [tstart tstop]
-    ## or [t0, t1, ..., tn], y is the initial value vector "init"
-    counter = 1;
-    fig = figure ();
+  if (isempty (flag))
+    ## Default case, plot and return a value
+    idx += 1;
+    told(idx,1) = t(1,1);
+    yold(:,idx) = y(:,1);
+    for i = 1:num_lines
+      set (hlines(i), "xdata", told, "ydata", yold(i,:));
+    endfor
+    drawnow;
+
+    retval = false;
+
+  elseif (strcmp (flag, "init"))
+    ## t is either the time slot [tstart tstop] or [t0, t1, ..., tn]
+    ## y is the initial value vector for the ode solution
+    idx = 1;
     told = t(1,1);
     yold = y(:,1);
-
-  elseif (isempty (flag))
-    ## Return something, either false for "not stopping
-    ## the integration" or true for "stopping the integration"
-    counter += 1;
-    figure (fig);
-    told(counter,1) = t(1,1);
-    yold(:,counter) = y(:,1);
-    ## FIXME: Why not use '.' rather than 'o' and skip the markersize?
-    ## FIXME: Why not just update the xdata, ydata properties?
-    ##        Calling plot involves a lot of overhead.
-    plot (told, yold, "-o", "markersize", 1); drawnow;
-    ret = false;
+    figure ();
+    hlines = plot (told, yold, "-", "marker", ".", "markersize", 9);
+    xlim ([t(1), t(end)]);  # Fix limits which also speeds up plotting
+    num_lines = numel (hlines);
 
   elseif (strcmp (flag, "done"))
     ## Cleanup after ode solver has finished.
-    clear ("figure", "told", "yold", "counter");
+    hlines = num_lines = told = yold = idx = [];
 
   endif
 
@@ -100,9 +113,9 @@
 
 
 %!demo
-%! # Solve an anonymous implementation of the Van der Pol equation
-%! # and display the results while solving
-%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
+%! ## Solve an anonymous implementation of the Van der Pol equation
+%! ## and display the results while solving
+%! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
 %! opt = odeset ("OutputFcn", @odeplot, "RelTol", 1e-6);
-%! sol = ode45 (fvdb, [0 20], [2 0], opt);
+%! sol = ode45 (fvdp, [0 20], [2 0], opt);
 
--- a/scripts/ode/odeset.m	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/odeset.m	Wed Oct 19 15:20:33 2016 -0700
@@ -87,10 +87,6 @@
 ## Mass matrix, specified as a constant matrix or a function of
 ## time and state.
 ##
-## @item MassConstant
-## Specify whether the mass matrix is a constant matrix or depends on
-## the state.
-##
 ## @item MassSingular
 ## Specify whether the mass matrix is singular. Accepted values include
 ## @qcode{"yes"}, @qcode{"no"}, @qcode{"maybe"}.
@@ -165,8 +161,6 @@
     p.addParameter ("JConstant", []);
     p.addParameter ("JPattern", []);
     p.addParameter ("Mass", []);
-    ## FIXME: MassConstant does not appear in Matlab documentation for odeset
-    p.addParameter ("MassConstant", []);
     p.addParameter ("MassSingular", []);
     p.addParameter ("MaxOrder", []);
     p.addParameter ("MaxStep", []);
@@ -190,22 +184,17 @@
     odestruct = p.Results;
     odestruct_extra = p.Unmatched;
 
-    ## FIXME: For speed, shouldn't this merge of structures only occur
-    ##        when there is something in odestruct_extra?
-    ## FIXME: Should alphabetical order of fieldnames be maintained
-    ##        by using sort?
-    s1 = cellfun (@(x) ifelse (iscell (x), {x}, x),
-                  struct2cell (odestruct),
-                  "UniformOutput", false);
+    xtra_fields = fieldnames (odestruct_extra);
+    if (! isempty (xtra_fields))
+      ## Merge extra fields into existing odestruct
+      for fldname = sort (xtra_fields.')
+        fldname = fldname{1};
+        warning ("Octave:invalid-input-arg",
+                 "odeset: unknown option \"%s\"\n", fldname);
+        odestruct.(fldname) = odestruct_extra.(fldname);
+      endfor
+    endif
 
-    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
@@ -225,7 +214,6 @@
   disp ('          JConstant:  binary, {["off"], "on"}');
   disp ('           JPattern:  sparse matrix, []');
   disp ('               Mass:  matrix or function_handle, []');
-  disp ('       MassConstant:  binary, {["off"], "on"}');
   disp ('       MassSingular:  switch, {["maybe"], "no", "yes"}');
   disp ('           MaxOrder:  switch, {[5], 1, 2, 3, 4, }');
   disp ('            MaxStep:  scalar, >0, []');
@@ -244,19 +232,19 @@
 
 
 %!demo
-%! # A new ODE options structure with default values is created.
+%! ## A new ODE options structure with default values is created.
 %!
 %! odeoptA = odeset ();
 
 %!demo
-%! # A new ODE options structure with manually set options
-%! # for "AbsTol" and "RelTol" is created.
+%! ## A new ODE options structure with manually set options
+%! ## for "AbsTol" and "RelTol" is created.
 %!
 %! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
 
 %!demo
-%! # A new ODE options structure is created from odeoptB with
-%! # a modified value for option "NormControl".
+%! ## A new ODE options structure is created from odeoptB with
+%! ## a modified value for option "NormControl".
 %!
 %! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1);
 %! odeoptC = odeset (odeoptB, "NormControl", "on");
@@ -264,7 +252,7 @@
 %!test
 %! odeoptA = odeset ();
 %! assert (isstruct (odeoptA));
-%! assert (numfields (odeoptA), 23);
+%! assert (numfields (odeoptA), 22);
 %! assert (all (structfun ("isempty", odeoptA)));
 
 %!shared odeoptB, odeoptC
@@ -284,13 +272,9 @@
 
 ## Test custom user-defined option
 %!test
-%! wstate = warning ("off", "Octave:invalid-input-arg");
-%! unwind_protect
-%!   odeopt = odeset ("NewtonTol", 3);
-%!   assert (odeopt.NewtonTol, 3);
-%! unwind_protect_cleanup
-%!   warning (wstate);
-%! end_unwind_protect
+%! warning ("off", "Octave:invalid-input-arg", "local");
+%! odeopt = odeset ("NewtonTol", 3);
+%! assert (odeopt.NewtonTol, 3);
 
 ## FIXME: Add an inexact match option once it is available in inputParser.
 ## See bug #49364.
@@ -299,12 +283,7 @@
 
 ## Test input validation
 %!error <argument 'OPT1' is not a valid parameter> odeset ("opt1")
-%!error  odeset (1, 1)
+%!error odeset (1, 1)
 %!error <argument 'OPT1' is not a valid parameter> odeset (odeset (), "opt1")
-%!error  odeset (odeset (), 1, 1)
+%!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	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/private/integrate_adaptive.m	Wed Oct 19 15:20:33 2016 -0700
@@ -96,8 +96,8 @@
     else
       solution.retout = x;
     endif
-    feval (options.OutputFcn, tspan, solution.retout,
-           "init", options.funarguments{:});
+    feval (options.OutputFcn, tspan, solution.retout, "init",
+           options.funarguments{:});
   endif
 
   ## Initialize the EventFcn
@@ -116,6 +116,7 @@
   solution.unhandledtermination = true;
   ireject = 0;
 
+  NormControl = strcmp (options.NormControl, "on");
   k_vals = [];
   iout = istep = 1;
 
@@ -133,9 +134,8 @@
       x_est(nn, end) = abs (x_est(nn, end));
     endif
 
-    ## FIXME: Take strcmp out of while loop and calculate just once
     err = AbsRel_norm (x_new, x_old, options.AbsTol, options.RelTol,
-                       strcmp (options.NormControl, "on"), x_est);
+                       NormControl, x_est);
 
     ## Accept solution only if err <= 1.0
     if (err <= 1)
@@ -186,7 +186,7 @@
           endif
 
           ## Call OutputFcn only if a valid result has been found.
-          ## Stop integration if function returns false.
+          ## Stop integration if function returns true.
           if (options.haveoutputfunction)
             cnt = options.Refine + 1;
             approxtime = linspace (t_old, t_new, cnt);
@@ -196,12 +196,16 @@
             if (! isempty (options.OutputSel))
               approxvals = approxvals(options.OutputSel, :);
             endif
+            stop_solve = false;
             for ii = 1:numel (approxtime)
-              pltret = feval (options.OutputFcn, approxtime(ii),
-                              approxvals(:, ii), [],
-                              options.funarguments{:});
+              stop_solve = feval (options.OutputFcn,
+                                  approxtime(ii), approxvals(:, ii), [],
+                                  options.funarguments{:});
+              if (stop_solve)
+                break;  # break from inner loop
+              endif
             endfor
-            if (pltret)  # Leave main loop
+            if (stop_solve)  # Leave main loop
               solution.unhandledtermination = false;
               break;
             endif
@@ -230,7 +234,7 @@
         endif
 
         ## Call OutputFcn only if a valid result has been found.
-        ## Stop integration if function returns false.
+        ## Stop integration if function returns true.
         if (options.haveoutputfunction)
           cnt = options.Refine + 1;
           approxtime = linspace (t_old, t_new, cnt);
@@ -240,11 +244,16 @@
           if (! isempty (options.OutputSel))
             approxvals = approxvals(options.OutputSel, :);
           endif
+          stop_solve = false;
           for ii = 1:numel (approxtime)
-            pltret = feval (options.OutputFcn, approxtime(ii),
-                            approxvals(:, ii), [], options.funarguments{:});
+            stop_solve = feval (options.OutputFcn,
+                                approxtime(ii), approxvals(:, ii), [],
+                                options.funarguments{:});
+            if (stop_solve)
+              break;  # break from inner loop
+            endif
           endfor
-          if (pltret)  # Leave main loop
+          if (stop_solve)  # Leave main loop
             solution.unhandledtermination = false;
             break;
           endif
@@ -280,7 +289,7 @@
     err += eps;  # avoid divisions by zero
     dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1))));
     dt = dir * min (abs (dt), options.MaxStep);
-    if (! (abs (dt) > eps (t (end))))
+    if (! (abs (dt) > eps (t(end))))
       break;
     endif
 
--- a/scripts/ode/private/ode_event_handler.m	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/private/ode_event_handler.m	Wed Oct 19 15:20:33 2016 -0700
@@ -19,44 +19,44 @@
 ## -*- texinfo -*-
 ## @deftypefn {} {@var{retval} =} ode_event_handler (@var{@@evtfun}, @var{t}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{})
 ##
-## Return the solution of the event function that is specified as the first
-## input argument @var{@@evtfun} in the form of a function handle.
+## Return the solution of the event function (@var{@@evtfun}) which is
+## specified in the form of a function handle.
+#
+## The second input argument @var{t} is a scalar double and specifies the time
+## of the event evaluation.
 ##
-## The second input argument @var{t} is of type double scalar and
-## specifies the time of the event evaluation, the third input argument
-## @var{y} either is of type double column vector (for ODEs and DAEs) and
-## specifies the solutions or is of type cell array (for IDEs and DDEs) and
-## specifies the derivatives or the history values, the third input argument
-## @var{flag} is of type string and can be of the form
+## The third input argument @var{y} may be a column vector of type double
+## (for ODEs and DAEs) which specifies the solutions.  Alternatives, @var{y}
+## may be a cell array (for IDEs and DDEs) which specifies the solutions and
+## derivatives.
+##
+## The fourth input argument @var{flag} is of type string.  Valid values are:
 ##
 ## @table @option
 ## @item  @qcode{"init"}
-## then initialize internal persistent variables of the function
-## @code{ode_event_handler} and return an empty cell array of size 4,
+## Initialize internal persistent variables of the function
+## @code{ode_event_handler} and return an empty cell array of size 4.
 ##
 ## @item  @qcode{"calc"}
-## then do the evaluation of the event function and return the solution
-## @var{retval} as type cell array of size 4,
+## Evaluate the event function and return the solution @var{retval} as a cell
+## array of size 4.
 ##
 ## @item  @qcode{"done"}
-## then cleanup internal variables of the function
-## @code{ode_event_handler} and return an empty cell array of size 4.
+## Clean up internal variables of the function @code{ode_event_handler} and
+## return an empty cell array of size 4.
 ## @end table
 ##
-## Optionally if further input arguments @var{par1}, @var{par2}, @dots{} of
-## any type are given then pass these parameters through
-## @code{ode_event_handler} to the event function.
+## If additional input arguments @var{par1}, @var{par2}, @dots{} are given
+## these parameters are passed directly to the event function.
 ##
-## This function is an ODE internal helper function therefore it should
-## never be necessary that this function is called directly by a user.  There
-## is only little error detection implemented in this function file to
-## achieve the highest performance.
+## This function is an ODE internal helper function and it should never be
+## necessary to call it directly.
 ## @end deftypefn
 
 function retval = ode_event_handler (evtfun, t, y, flag = "", varargin)
 
   ## No error handling has been implemented in this function to achieve
-  ## the highest performance available.
+  ## the highest performance possible.
 
   ## retval{1} is true (to terminate) or false (to continue)
   ## retval{2} is the index information for which an event occurred
@@ -70,33 +70,12 @@
   ## yold    the ODE result
   ## retcell the return values cell array
   ## evtcnt  the counter for how often this function has been called
-  persistent evtold told yold retcell evtcnt;
-
-  ## Call the event function if an event function has been defined to
-  ## initialize the internal variables of the event function and to get
-  ## a value for evtold.
-  if (strcmp (flag, "init"))
+  persistent evtold told yold retcell;
+  persistent evtcnt = 1;   # Don't remove.  Required for Octave parser.
 
-    if (! iscell (y))
-      inpargs = {evtfun, t, y};
-    else
-      inpargs = {evtfun, t, y{1}, y{2}};
-      y = y{1};  # Delete cell element 2
-    endif
-    if (nargin > 4)
-      inpargs = {inpargs{:}, varargin{:}};
-    endif
-    [evtold, term, dir] = feval (inpargs{:});
-
-    ## FIXME: This actually seems to assume that everything must be row vectors
-    ## We assume that all return values must be column vectors
-    evtold = evtold(:)'; term = term(:)'; dir = dir(:)';
-    told = t; yold = y; evtcnt = 1; retcell = cell (1,4);
-
-  ## Process the event, i.e.,
-  ## find the zero crossings for either a rising or falling edge
-  elseif (isempty (flag))
-
+  if (isempty (flag))
+    ## Process the event, i.e.,
+    ## find the zero crossings for either a rising or falling edge
     if (! iscell (y))
       inpargs = {evtfun, t, y};
     else
@@ -108,8 +87,7 @@
     endif
     [evt, term, dir] = feval (inpargs{:});
 
-    ## FIXME: This actually seems to assume that everything must be row vectors
-    ## We assume that all return values must be column vectors
+    ## We require that all return values be row vectors
     evt = evt(:)'; term = term(:)'; dir = dir(:)';
 
     ## Check if one or more signs of the event has changed
@@ -143,12 +121,35 @@
         evtcnt += 1;
       endif
 
-    endif  # Check for one or more signs ...
-    evtold = evt; told = t; retval = retcell; yold = y;
+    endif
+    evtold = evt; told = t; yold = y;
+    retval = retcell;
+
+  elseif (strcmp (flag, "init"))
+    ## Call the event function if an event function has been defined to
+    ## initialize the internal variables of the event function and to get
+    ## a value for evtold.
 
-  elseif (strcmp (flag, "done"))  # Clear this event handling function
-    clear ("evtold", "told", "retcell", "yold", "evtcnt");
-    retcell = cell (1,4);
+    if (! iscell (y))
+      inpargs = {evtfun, t, y};
+    else
+      inpargs = {evtfun, t, y{1}, y{2}};
+      y = y{1};  # Delete cell element 2
+    endif
+    if (nargin > 4)
+      inpargs = {inpargs{:}, varargin{:}};
+    endif
+    [evtold, ~, ~] = feval (inpargs{:});
+
+    ## We require that all return values be row vectors
+    evtold = evtold(:)'; told = t; yold = y;
+    evtcnt = 1;
+    retval = retcell = cell (1,4);
+
+  elseif (strcmp (flag, "done"))
+    ## Clear this event handling function
+    evtold = told = yold = evtcnt = [];
+    retval = retcell = cell (1,4);
 
   endif
 
--- a/scripts/ode/private/odedefaults.m	Mon Oct 17 16:26:10 2016 -0400
+++ b/scripts/ode/private/odedefaults.m	Wed Oct 19 15:20:33 2016 -0700
@@ -33,10 +33,9 @@
                                 "JConstant", "off",
                                 "JPattern", [],
                                 "Mass", [],
-                                "MassConstant", "off",
                                 "MassSingular", "maybe",
                                 "MaxOrder", 5,
-                                "MaxStep", 0.1 * abs (t0-tf),
+                                "MaxStep", 0.1 * abs (tf - t0),
                                 "MStateDependence", "weak",
                                 "MvPattern", [],
                                 "NonNegative", [],
@@ -48,7 +47,7 @@
                                 "Stats", "off",
                                 "Vectorized", "off");
 
-  defaults.MaxStep = (0.1 * abs (t0-tf));
+  defaults.MaxStep = 0.1 * abs (tf -t0);
 
   persistent classes = struct ("AbsTol", {{"float"}},
                                "BDF", "char",
@@ -59,7 +58,6 @@
                                "JConstant", "char",
                                "JPattern", {{"float"}},
                                "Mass", {{"float", "function_handle"}},
-                               "MassConstant", "char",
                                "MassSingular", "char",
                                "MaxOrder", {{"float"}},
                                "MaxStep", {{"float"}},
@@ -83,7 +81,6 @@
                                   "JConstant", {{"on", "off"}},
                                   "JPattern", {{"vector"}},
                                   "Mass", {{}},
-                                  "MassConstant", {{"on", "off"}},
                                   "MassSingular", {{"no", "maybe", "yes"}},
                                   "MaxOrder", {{">=", 0, "<=", 5, "integer"}},
                                   "MaxStep", {{"positive", "scalar", "real"}},