# HG changeset patch
# User Ken Marek <marek_ka@mercer.edu>
# Date 1664322516 14400
#      Tue Sep 27 19:48:36 2022 -0400
# Node ID 2bd7713ccdb8a0c5e07a14d0ceebe7d93681ab3a
# Parent  dd6b37f67db2939f39cd99717a69b64e917cbc7e
ode45, ode23, ode23s: Events, OutputFcn, & Refine options improvement/implementation (bug #63063) (subsumes bug #49408). Improves Matlab compatibility with results of these functions.

* scripts/ode/ode23.m: Various changes for Matlab compatibility and due to structural changes in integrate_adaptive.m. Remove disabling of Refine option when outputting struct (Matlab compatibility: still use Refine for OutputFcn calls). Transpose solution variables when outputting struct (Matlab compatibility). Solution struct now has two sets of solution variables: output_t, output_x and ode_t and ode_x; this reflects new data handling in integrate_adaptive.m for Matlab compatibility. Update file tests to reflect changed functionality. Minor code formatting.
* scripts/ode/ode23s.m: Same changes as ode23.m.
* scripts/ode/ode45.m: Same changes as ode23.m. Also remove comment indicating that Refine is not implemented.
* scripts/ode/private/integrate_adaptive.m: Handing of Refine, Events, OutputFcn, and time output with specified times are all changed to improve Matlab compatibility. Update internal handling of variables t and x, including keeping separate records for internal integration use (ode_t, ode_x) and potentially interpolated results for function output or calls to OutputFcn by this script (output_t, output_x). Apparent old attempt at implementing Refine option is removed and replaced with new implementation. Specified time output or Refine != 0 are both interpolated from internal variables (ode_t, ode_x) for output of non-struct variables and/or for use with OutputFcn. Event handling is improved where multiple Events (including at least one terminal Event) are detected in a single simulation step, so that all Events (up to and including the first terminal one) are reported, and final data point is set to that of terminal Event. Multiple data points are sent in a single call to OutputFcn if they are all interpolated from a single integration step (Matlab compatibility). Remove warning for termination when termination signal is received from Events or OutputFcn (Matlab compatibility). Return internal variables (ode_t, ode_x) and interpolated variables (output_t, output_x) so calling function can correctly return struct or separate variables as needed (Matlab compatibility).
* scripts/ode/private/ode_event_handler.m: Handing of Events is improved for better Matlab compatibility. When multiple Events are encountered in one integration step, they are sorted in ascending time order. If there is at least one terminal Event, all Events after the time of the terminal Event are removed. These changes also allow integrate_adaptive.m to assume that when a terminal Event is recorded, the last Event that is returned is at the terminal time.
* scripts/ode/odeset.m: Update docstring to remove indication that Refine is not implemented
* scripts/ode/odeplot.m: Update docstring to indicate that input t can be a scalar or vector. The code seemed to work properly with multiple time inputs in one function call, so no changes were made to functional code. Added file test. This file needed to be checked/updated because it is the default function for OutputFcn, and needs to be compatible with updated OutputFcn option behavior.
* etc/NEWS.8.md: Add descriptions of changes under General improvements and Matlab compatibility.

diff -r dd6b37f67db2 -r 2bd7713ccdb8 etc/NEWS.8.md
--- a/etc/NEWS.8.md	Sun Sep 25 06:22:25 2022 -0400
+++ b/etc/NEWS.8.md	Tue Sep 27 19:48:36 2022 -0400
@@ -22,6 +22,9 @@
 - `quadgk` can now accept the `ArrayValued` input parameter to handle
 array-valued input functions.
 
+- `Refine` option is now implemented in functions `ode45`, `ode23`, 
+  and `ode23s`.
+
 ### Graphical User Interface
 
 
@@ -68,6 +71,10 @@
   Object      | Property         | Default State
   ------------|------------------|------------
   `figure`    | `"dockcontrols"` | `"on"`
+  
+- `ode45`, `ode23`, and `ode23s` have improved results for options `Events`,
+  `OutputFcn`, and `Refine`, along with corrected orientation of struct 
+  outputs.
 
 ### Alphabetical list of new functions added in Octave 8
 
diff -r dd6b37f67db2 -r 2bd7713ccdb8 scripts/ode/ode23.m
--- a/scripts/ode/ode23.m	Sun Sep 25 06:22:25 2022 -0400
+++ b/scripts/ode/ode23.m	Tue Sep 27 19:48:36 2022 -0400
@@ -196,8 +196,8 @@
     odeopts.InitialStep = odeopts.direction * ...
                           starting_stepsize (order, fcn, trange(1), init,
                                              odeopts.AbsTol, odeopts.RelTol,
-                                             strcmpi (odeopts.NormControl, "on"),
-                                             odeopts.funarguments);
+                                             strcmpi (odeopts.NormControl,
+                                             "on"), odeopts.funarguments);
   endif
 
   if (! isempty (odeopts.Mass))
@@ -229,12 +229,7 @@
     endif
   endif
 
-  if (nargout == 1)
-    ## Single output requires auto-selected intermediate times,
-    ## which is obtained by NOT specifying specific solution times.
-    trange = [trange(1); trange(end)];
-    odeopts.Refine = [];  # disable Refine when single output requested
-  elseif (numel (trange) > 2)
+  if (numel (trange) > 2)
     odeopts.Refine = [];  # disable Refine when specific times requested
   endif
 
@@ -246,8 +241,9 @@
     feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
-    ode_event_handler (odeopts.Events, solution.t(end),
-                       solution.x(end,:).', "done", odeopts.funarguments{:});
+    ode_event_handler (odeopts.Events, solution.ode_t(end), ...
+                       solution.ode_x(end,:).', "done", ...
+                       odeopts.funarguments{:});
   endif
 
   ## Print additional information if option Stats is set
@@ -265,16 +261,16 @@
   endif
 
   if (nargout == 2)
-    varargout{1} = solution.t;      # Time stamps are first output argument
-    varargout{2} = solution.x;      # Results are second output argument
+    varargout{1} = solution.output_t; # Time stamps are first output argument
+    varargout{2} = solution.output_x; # Results are second output argument
   elseif (nargout == 1)
-    varargout{1}.x = solution.t.';  # Time stamps saved in field x (row vector)
-    varargout{1}.y = solution.x.';  # Results are saved in field y (row vector)
-    varargout{1}.solver = solver;   # Solver name is saved in field solver
+    varargout{1}.x = solution.ode_t.'; #Time stamps saved in field x (row vect.)
+    varargout{1}.y = solution.ode_x.'; #Results are saved in field y (row vect.)
+    varargout{1}.solver = solver; # Solver name is saved in field solver
     if (! isempty (odeopts.Events))
-      varargout{1}.xe = solution.event{3};  # Time info when an event occurred
-      varargout{1}.ye = solution.event{4};  # Results when an event occurred
-      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
+      varargout{1}.ie = solution.event{2}.'; # Index info which event occurred
     endif
     if (strcmpi (odeopts.Stats, "on"))
       varargout{1}.stats = struct ();
@@ -287,8 +283,8 @@
     endif
   elseif (nargout > 2)
     varargout = cell (1,5);
-    varargout{1} = solution.t;
-    varargout{2} = solution.x;
+    varargout{1} = solution.output_t;
+    varargout{2} = solution.output_x;
     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
@@ -355,7 +351,8 @@
 %!      error ('fout: step "init"');
 %!    endif
 %!  elseif (isempty (flag))
-%!    if (! isequal (size (t), [1, 1]))
+%!  # Multiple steps can be sent in one function call
+%!    if (! isequal ( size (t), size (y)))
 %!      error ('fout: step "calc"');
 %!    endif
 %!  elseif (strcmp (flag, "done"))
@@ -370,6 +367,14 @@
 %!test  # two output arguments
 %! [t, y] = ode23 (@fpol, [0 2], [2 0]);
 %! assert ([t(end), y(end,:)], [2, fref], 1e-3);
+%!test  # correct number of steps with Refine
+%! [t1, y1] = ode23 (@fpol, [0 2], [2 0], odeset ("Refine", 1));
+%! [t2, y2] = ode23 (@fpol, [0 2], [2 0], odeset ("Refine", 4));
+%! [t3, y3] = ode23 (@fpol, [0 2], [2 0]); #default Refine=1
+%! s = ode23 (@fpol, [0 2], [2 0], odeset ("Refine", 4));
+%! assert (length (t1) == length (t3));
+%! assert (length (t2) == 4*length (t1) - 3);
+%! assert (length (s.x) == length (t1));
 %!test  # anonymous function instead of real function
 %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
 %! [t, y] = ode23 (fvdp, [0 2], [2 0]);
@@ -435,8 +440,8 @@
 %! opt = odeset ("NonNegative", 2);
 %! sol = ode23 (@fpol, [0 2], [2 0], opt);
 %! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 1e-1);
-%!test  # Details of OutputSel and Refine can't be tested
-%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
+%!test  # Details of OutputSel can't be tested
+%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1);
 %! sol = ode23 (@fpol, [0 2], [2 0], opt);
 %!test  # Stats must add further elements in sol
 %! opt = odeset ("Stats", "on");
@@ -452,13 +457,11 @@
 %! assert (isfield (sol, "xe"));
 %! assert (isfield (sol, "ye"));
 %!test  # Events option, now stop integration
-%! warning ("off", "integrate_adaptive:unexpected_termination", "local");
 %! opt = odeset ("Events", @fevn, "NormControl", "on");
 %! sol = ode23 (@fpol, [0 10], [2 0], opt);
-%! assert ([sol.ie, sol.xe, sol.ye],
+%! assert ([sol.ie, sol.xe, sol.ye.'],
 %!         [2.0, 2.496110, -0.830550, -2.677589], .5e-1);
 %!test  # Events option, five output arguments
-%! warning ("off", "integrate_adaptive:unexpected_termination", "local");
 %! opt = odeset ("Events", @fevn, "NormControl", "on");
 %! [t, y, vxe, ye, vie] = ode23 (@fpol, [0 10], [2 0], opt);
 %! assert ([vie, vxe, ye], [2.0, 2.496110, -0.830550, -2.677589], 1e-1);
@@ -500,6 +503,9 @@
 %! assert (imag (sol.y), ones (size (sol.y)));
 %! [x, y] = ode23 (@(x,y) 1, [0 1], 1i);
 %! assert (imag (y), ones (size (y)));
+%!test # Make sure odeplot works (default OutputFcn when no return value)
+%! ode23 (@fpol, [0 2], [2 0]);
+%! close all
 
 ## Test input validation
 %!error <Invalid call> ode23 ()
diff -r dd6b37f67db2 -r 2bd7713ccdb8 scripts/ode/ode23s.m
--- a/scripts/ode/ode23s.m	Sun Sep 25 06:22:25 2022 -0400
+++ b/scripts/ode/ode23s.m	Tue Sep 27 19:48:36 2022 -0400
@@ -205,12 +205,7 @@
 
   ## Starting the initialization of the core solver ode23s
 
-  if (nargout == 1)
-    ## Single output requires auto-selected intermediate times,
-    ## which is obtained by NOT specifying specific solution times.
-    trange = [trange(1); trange(end)];
-    odeopts.Refine = [];  # disable Refine when single output requested
-  elseif (numel (trange) > 2)
+  if (numel (trange) > 2)
     odeopts.Refine = [];  # disable Refine when specific times requested
   endif
 
@@ -222,8 +217,9 @@
     feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
-    ode_event_handler (odeopts.Events, solution.t(end), ...
-                       solution.x(end,:).', "done", odeopts.funarguments{:});
+    ode_event_handler (odeopts.Events, solution.ode_t(end), ...
+                       solution.ode_x(end,:).', "done", ...
+                       odeopts.funarguments{:});
   endif
 
   ## Print additional information if option Stats is set
@@ -241,16 +237,16 @@
   endif
 
   if (nargout == 2)
-    varargout{1} = solution.t;      # Time stamps are first output argument
-    varargout{2} = solution.x;      # Results are second output argument
+    varargout{1} = solution.output_t; # Time stamps are first output argument
+    varargout{2} = solution.output_x; # Results are second output argument
   elseif (nargout == 1)
-    varargout{1}.x = solution.t.';  # Time stamps saved in field x (row vector)
-    varargout{1}.y = solution.x.';  # Results are saved in field y (row vector)
+    varargout{1}.x = solution.ode_t.'; #Time stamps saved in field x (row vect.)
+    varargout{1}.y = solution.ode_x.'; #Results are saved in field y (row vect.)
     varargout{1}.solver = solver;   # Solver name is saved in field solver
     if (! isempty (odeopts.Events))
-      varargout{1}.xe = solution.event{3};  # Time info when an event occurred
-      varargout{1}.ye = solution.event{4};  # Results when an event occurred
-      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
+      varargout{1}.ie = solution.event{2}.'; # Index info which event occurred
     endif
     if (strcmpi (odeopts.Stats, "on"))
       varargout{1}.stats = struct ();
@@ -263,8 +259,8 @@
     endif
   elseif (nargout == 5)
     varargout = cell (1,5);
-    varargout{1} = solution.t;
-    varargout{2} = solution.x;
+    varargout{1} = solution.output_t;
+    varargout{2} = solution.output_x;
     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
@@ -389,7 +385,8 @@
 %!      error ('fout: step "init"');
 %!    endif
 %!  elseif (isempty (flag))
-%!    if (! isequal (size (t), [1, 1]))
+%!  # Multiple steps can be sent in one function call
+%!    if (! isequal ( size (t), size (y)))
 %!      error ('fout: step "calc"');
 %!    endif
 %!  elseif (strcmp (flag, "done"))
@@ -404,6 +401,14 @@
 %!test  # two output arguments
 %! [t, y] = ode23s (@fpol, [0 2], [2 0]);
 %! assert ([t(end), y(end,:)], [2, fref], 1e-3);
+%!test  # correct number of steps with Refine
+%! [t1, y1] = ode23s (@fpol, [0 2], [2 0], odeset ("Refine", 1));
+%! [t2, y2] = ode23s (@fpol, [0 2], [2 0], odeset ("Refine", 4));
+%! [t3, y3] = ode23s (@fpol, [0 2], [2 0]); #default Refine=1
+%! s = ode23s (@fpol, [0 2], [2 0], odeset ("Refine", 4));
+%! assert (length (t1) == length (t3));
+%! assert (length (t2) == 4*length (t1) - 3);
+%! assert (length (s.x) == length (t1));
 %!test  # anonymous function instead of real function
 %! fvdp = @(t,y) [y(2); 10 * (1 - y(1)^2) * y(2) - y(1)];
 %! [t, y] = ode23s (fvdp, [0 2], [2 0]);
@@ -435,8 +440,8 @@
 %! opt = odeset ("RelTol", 1e-8, "NormControl", "on");
 %! sol = ode23s (@fpol, [0 2], [2 0], opt);
 %! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-4);
-%!test  # Details of OutputSel and Refine can't be tested
-%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
+%!test  # Details of OutputSel can't be tested
+%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1);
 %! sol = ode23s (@fpol, [0 2], [2 0], opt);
 %!test  # Stats must add further elements in sol
 %! opt = odeset ("Stats", "on");
@@ -498,6 +503,9 @@
 %! assert (imag (sol.y), ones (size (sol.y)));
 %! [x, y] = ode23s (@(x,y) 1, [0 1], 1i);
 %! assert (imag (y), ones (size (y)));
+%!test # Make sure odeplot works (default OutputFcn when no return value)
+%! ode23s (@fpol, [0 2], [2 0]);
+%! close all
 
 ## Test input validation
 %!error <Invalid call> ode23s ()
diff -r dd6b37f67db2 -r 2bd7713ccdb8 scripts/ode/ode45.m
--- a/scripts/ode/ode45.m	Sun Sep 25 06:22:25 2022 -0400
+++ b/scripts/ode/ode45.m	Tue Sep 27 19:48:36 2022 -0400
@@ -156,7 +156,6 @@
   [defaults, classes, attributes] = odedefaults (numel (init),
                                                  trange(1), trange(end));
 
-  ## FIXME: Refine is not correctly implemented yet
   defaults = odeset (defaults, "Refine", 4);
 
   persistent ode45_ignore_options = ...
@@ -196,8 +195,8 @@
     odeopts.InitialStep = odeopts.direction * ...
                           starting_stepsize (order, fcn, trange(1), init,
                                              odeopts.AbsTol, odeopts.RelTol,
-                                             strcmpi (odeopts.NormControl, "on"),
-                                             odeopts.funarguments);
+                                             strcmpi (odeopts.NormControl,
+                                             "on"), odeopts.funarguments);
   endif
 
   if (! isempty (odeopts.Mass))
@@ -229,12 +228,7 @@
     endif
   endif
 
-  if (nargout == 1)
-    ## Single output requires auto-selected intermediate times,
-    ## which is obtained by NOT specifying specific solution times.
-    trange = [trange(1); trange(end)];
-    odeopts.Refine = [];  # disable Refine when single output requested
-  elseif (numel (trange) > 2)
+  if (numel (trange) > 2)
     odeopts.Refine = [];  # disable Refine when specific times requested
   endif
 
@@ -246,8 +240,9 @@
     feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
-    ode_event_handler (odeopts.Events, solution.t(end),
-                       solution.x(end,:).', "done", odeopts.funarguments{:});
+    ode_event_handler (odeopts.Events, solution.ode_t(end), ...
+                       solution.ode_x(end,:).', "done", ...
+                       odeopts.funarguments{:});
   endif
 
   ## Print additional information if option Stats is set
@@ -265,16 +260,16 @@
   endif
 
   if (nargout == 2)
-    varargout{1} = solution.t;      # Time stamps are first output argument
-    varargout{2} = solution.x;      # Results are second output argument
+    varargout{1} = solution.output_t; # Time stamps are first output argument
+    varargout{2} = solution.output_x; # Results are second output argument
   elseif (nargout == 1)
-    varargout{1}.x = solution.t.';  # Time stamps saved in field x (row vector)
-    varargout{1}.y = solution.x.';  # Results are saved in field y (row vector)
-    varargout{1}.solver = solver;   # Solver name is saved in field solver
+    varargout{1}.x = solution.ode_t.'; #Time stamps saved in field x (row vect.)
+    varargout{1}.y = solution.ode_x.'; #Results are saved in field y (row vect.)
+    varargout{1}.solver = solver; # Solver name is saved in field solver
     if (! isempty (odeopts.Events))
-      varargout{1}.xe = solution.event{3};  # Time info when an event occurred
-      varargout{1}.ye = solution.event{4};  # Results when an event occurred
-      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
+      varargout{1}.ie = solution.event{2}.'; # Index info which event occurred
     endif
     if (strcmpi (odeopts.Stats, "on"))
       varargout{1}.stats = struct ();
@@ -287,8 +282,8 @@
     endif
   elseif (nargout > 2)
     varargout = cell (1,5);
-    varargout{1} = solution.t;
-    varargout{2} = solution.x;
+    varargout{1} = solution.output_t;
+    varargout{2} = solution.output_x;
     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
@@ -355,7 +350,8 @@
 %!      error ('fout: step "init"');
 %!    endif
 %!  elseif (isempty (flag))
-%!    if (! isequal (size (t), [1, 1]))
+%!  # Multiple steps can be sent in one function call
+%!    if (! isequal ( size (t), size (y)))
 %!      error ('fout: step "calc"');
 %!    endif
 %!  elseif (strcmp (flag, "done"))
@@ -371,8 +367,16 @@
 %! [t, y] = ode45 (@fpol, [0 2], [2 0]);
 %! assert ([t(end), y(end,:)], [2, fref], 1e-2);
 %!test  # not too many steps
-%! [t, y] = ode45 (@fpol, [0 2], [2 0]);
+%! [t, y] = ode45 (@fpol, [0 2], [2 0], odeset("Refine",1));
 %! assert (size (t) < 20);
+%!test  # correct number of steps with Refine
+%! [t1, y1] = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 1));
+%! [t2, y2] = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 4));
+%! [t3, y3] = ode45 (@fpol, [0 2], [2 0]); #default Refine=4
+%! s = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 4));
+%! assert (length (t2) == length (t3));
+%! assert (length (t2) == 4*length (t1) - 3);
+%! assert (length (s.x) == length (t1));
 %!test  # anonymous function instead of real function
 %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
 %! [t, y] = ode45 (fvdp, [0 2], [2 0]);
@@ -392,7 +396,7 @@
 %! [t, y] = ode45 (@(t,y) y, [-2 0], 2);
 %! assert ([t(end), y(end,:)], vref, 1e-1);
 %!test  # InitialStep option
-%! opt = odeset ("InitialStep", 1e-8);
+%! opt = odeset ("InitialStep", 1e-8, "Refine", 1);
 %! [t, y] = ode45 (@fpol, [0 0.2], [2 0], opt);
 %! assert ([t(2)-t(1)], [1e-8], 1e-9);
 %!test  # MaxStep option
@@ -447,8 +451,8 @@
 %! opt = odeset ("NonNegative", 2);
 %! sol = ode45 (@fpol, [0 2], [2 0], opt);
 %! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 0.5);
-%!test  # Details of OutputSel and Refine can't be tested
-%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
+%!test  # Details of OutputSel can't be tested
+%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1);
 %! sol = ode45 (@fpol, [0 2], [2 0], opt);
 %!test  # Stats must add further elements in sol
 %! opt = odeset ("Stats", "on");
@@ -464,13 +468,11 @@
 %! assert (isfield (sol, "xe"));
 %! assert (isfield (sol, "ye"));
 %!test  # Events option, now stop integration
-%! warning ("off", "integrate_adaptive:unexpected_termination", "local");
 %! opt = odeset ("Events", @fevn, "NormControl", "on");
 %! sol = ode45 (@fpol, [0 10], [2 0], opt);
-%! assert ([sol.ie, sol.xe, sol.ye],
+%! assert ([sol.ie, sol.xe, sol.ye.'],
 %!         [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
 %!test  # Events option, five output arguments
-%! warning ("off", "integrate_adaptive:unexpected_termination", "local");
 %! opt = odeset ("Events", @fevn, "NormControl", "on");
 %! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt);
 %! assert ([vie, vxe, ye],
@@ -511,8 +513,11 @@
 %!test # Check that imaginary part of solution does not get inverted
 %! sol = ode45 (@(x,y) 1, [0 1], 1i);
 %! assert (imag (sol.y), ones (size (sol.y)));
-%! [x, y] = ode45 (@(x,y) 1, [0 1], 1i);
+%! [x, y] = ode45 (@(x,y) 1, [0 1], 1i, odeset ("Refine", 1));
 %! assert (imag (y), ones (size (y)));
+%!test # Make sure odeplot works (default OutputFcn when no return value)
+%! ode45 (@fpol, [0 2], [2 0]);
+%! close all
 
 %!error <Invalid call> ode45 ()
 %!error <Invalid call> ode45 (1)
diff -r dd6b37f67db2 -r 2bd7713ccdb8 scripts/ode/odeplot.m
--- a/scripts/ode/odeplot.m	Sun Sep 25 06:22:25 2022 -0400
+++ b/scripts/ode/odeplot.m	Tue Sep 27 19:48:36 2022 -0400
@@ -40,8 +40,8 @@
 ## contains the initial conditions for the ode problem (@var{y0}).
 ##
 ## @item @qcode{""}
-## The input @var{t} must be a scalar double specifying the time for which
-## the solution in input @var{y} was calculated.
+## The input @var{t} must be a scalar double or vector specifying the time(s)
+## for which the solution in input @var{y} was calculated.
 ##
 ## @item @qcode{"done"}
 ## The inputs should be empty, but are ignored if they are present.
@@ -120,3 +120,14 @@
 %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
 %! opt = odeset ("OutputFcn", @odeplot, "RelTol", 1e-6);
 %! sol = ode45 (fvdp, [0 20], [2 0], opt);
+
+## Test that the function works for expected ode OutputFcn calls
+
+%!test
+%! t = linspace(0,2,10);
+%! y = [0.2*t; -0.1*t.^2-1; sin(t)];
+%! odeplot ([0 2], y(:,1), "init");
+%! odeplot (t(2), y(:,2), []);
+%! odeplot (t(3:end), y(:, 3:end), []);
+%! odeplot ([], [], "done");
+%! close all
diff -r dd6b37f67db2 -r 2bd7713ccdb8 scripts/ode/odeset.m
--- a/scripts/ode/odeset.m	Sun Sep 25 06:22:25 2022 -0400
+++ b/scripts/ode/odeset.m	Tue Sep 27 19:48:36 2022 -0400
@@ -130,7 +130,6 @@
 ## time step or also at intermediate time instances.  The value should be
 ## a scalar indicating the number of equally spaced time points to use
 ## within each timestep at which to return output.
-## @emph{Note}: This option is not yet implemented.
 ##
 ## @item @code{RelTol}: positive scalar
 ## Relative error tolerance.
diff -r dd6b37f67db2 -r 2bd7713ccdb8 scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m	Sun Sep 25 06:22:25 2022 -0400
+++ b/scripts/ode/private/integrate_adaptive.m	Tue Sep 27 19:48:36 2022 -0400
@@ -73,13 +73,13 @@
 
   fixed_times = numel (tspan) > 2;
 
-  t_new = t_old = t = tspan(1);
-  x_new = x_old = x = x0(:);
+  t_new = t_old = ode_t = output_t = tspan(1);
+  x_new = x_old = ode_x = output_x = x0(:);
 
   ## Get first initial timestep
   dt = options.InitialStep;
   if (isempty (dt))
-    dt = starting_stepsize (order, fcn, t, x,
+    dt = starting_stepsize (order, fcn, ode_t, ode_x,
                             options.AbsTol, options.RelTol,
                             strcmp (options.NormControl, "on"),
                             options.funarguments);
@@ -95,12 +95,23 @@
   facmax = 1.5;
   fac = 0.38^(1/(order+1));  # formula taken from Hairer
 
+  ## Initialize Refine value
+  refine = options.Refine;
+  if isempty (refine)
+    refine = 1;
+  elseif ((refine != round (refine)) || (refine < 1))
+    refine = 1;
+    warning ("integrate_adaptive:invalid_refine",
+               ["Invalid value of Refine. Refine must be a positive " ...
+                "integer. Setting Refine = 1."] );
+  endif
+
   ## Initialize the OutputFcn
   if (options.haveoutputfunction)
     if (! isempty (options.OutputSel))
-      solution.retout = x(options.OutputSel,end);
+      solution.retout = output_x(options.OutputSel, end);
     else
-      solution.retout = x;
+      solution.retout = output_x;
     endif
     feval (options.OutputFcn, tspan, solution.retout, "init",
            options.funarguments{:});
@@ -110,7 +121,7 @@
   have_EventFcn = false;
   if (! isempty (options.Events))
     have_EventFcn  = true;
-    ode_event_handler (options.Events, tspan(1), x,
+    ode_event_handler (options.Events, tspan(1), ode_x,
                        "init", options.funarguments{:});
   endif
 
@@ -150,130 +161,84 @@
 
       solution.cntloop += 1;
       ireject = 0;              # Clear reject counter
-
-      ## if output time steps are fixed
-      if (fixed_times)
-
-        t_caught = find ((dir * tspan(iout:end) > dir * t_old)
-                         & (dir * tspan(iout:end) <= dir * t_new));
-        t_caught = t_caught + iout - 1;
+      terminal_event = false;
+      terminal_output = false;
 
-        if (! isempty (t_caught))
-          t(t_caught) = tspan(t_caught);
-          iout = max (t_caught);
-          x(:, t_caught) = ...
-            runge_kutta_interpolate (order, [t_old t_new], [x_old x_new],
-                                     t(t_caught), new_k_vals, dt, fcn,
-                                     options.funarguments);
-
-          istep += 1;
+      istep++;
+      ode_t(istep) = t_new;
+      ode_x(:, istep) = x_new;
+      iadd = 0;                 # Number of output points added this iteration
 
-          ## Call Events function only if a valid result has been found.
-          ## Stop integration if eventbreak is true.
-          if (have_EventFcn)
-            break_loop = false;
-            for idenseout = 1:numel (t_caught)
-              id = t_caught(idenseout);
-              td = t(id);
-              solution.event = ...
-                ode_event_handler (options.Events, t(id), x(:, id), [],
-                                   options.funarguments{:});
-              if (! isempty (solution.event{1}) && solution.event{1} == 1)
-                t(id) = solution.event{3}(end);
-                t = t(1:id);
-                x(:, id) = solution.event{4}(end, :).';
-                x = x(:,1:id);
-                solution.unhandledtermination = false;
-                break_loop = true;
-                break;
-              endif
-            endfor
-            if (break_loop)
-              break;
-            endif
-          endif
+      ## Check for Events
+      if (have_EventFcn)
+        solution.event = ode_event_handler (options.Events, t_new, x_new, ...
+                                            [], options.funarguments{:});
+        ## Check for terminal Event
+        if (! isempty (solution.event{1}) && solution.event{1} == 1)
+          ode_t(istep) = solution.event{3}(end);
+          ode_x(:, istep) = solution.event{4}(end, :).';
+          solution.unhandledtermination = false;
+          terminal_event = true;
+        endif
+      endif
 
-          ## Call OutputFcn only if a valid result has been found.
-          ## Stop integration if function returns true.
-          if (options.haveoutputfunction)
-            cnt = options.Refine + 1;
-            approxtime = linspace (t_old, t_new, cnt);
-            approxvals = interp1 ([t_old, t(t_caught), t_new],
-                                  [x_old, x(:, t_caught), x_new] .',
-                                  approxtime, "linear") .';
-            if (isvector (approxvals) && ! isscalar (approxtime))
-              approxvals = approxvals.';
-            endif
-            if (! isempty (options.OutputSel))
-              approxvals = approxvals(options.OutputSel, :);
-            endif
-            stop_solve = false;
-            for ii = 1:numel (approxtime)
-              stop_solve = feval (options.OutputFcn,
-                                  approxtime(ii), approxvals(:, ii), [],
-                                  options.funarguments{:});
-              if (stop_solve)
-                break;  # break from inner loop
-              endif
-            endfor
-            if (stop_solve)  # Leave main loop
-              solution.unhandledtermination = false;
-              break;
-            endif
-          endif
-
+      ## Interpolate to specified or Refined points
+      if (fixed_times)
+        t_caught = find ((dir * tspan(iout:end) > dir * t_old) ...
+                         & (dir * tspan(iout:end) <= dir * ode_t(istep)));
+        t_caught = t_caught + iout - 1;
+        iadd = length (t_caught);
+        if (! isempty (t_caught))
+          output_t(t_caught) = tspan(t_caught);
+          iout = max (t_caught);
+          output_x(:, t_caught) = ...
+            runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], ...
+                                     output_t(t_caught), new_k_vals, dt, ...
+                                     fcn, options.funarguments);
         endif
-
-      else   # not fixed times
-
-        t(++istep)  = t_new;
-        x(:, istep) = x_new;
-        iout = istep;
+        ## Add a possible additional output value if we found a terminal Event
+        if ((terminal_event == true) && ...
+            (dir * ode_t(istep) > dir * output_t(iout)))
+          iadd += 1;
+          iout += 1;
+          output_x(:, iout) = ode_x(:, istep);
+          output_t(iout) = ode_t(istep);
+        endif
+      elseif (refine > 1)
+        iadd = refine;
+        tadd = linspace (t_old, ode_t(istep), refine + 1);
+        tadd = tadd(2:end);
+        output_x(:, iout + (1:iadd)) = ...
+          runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], ...
+                                   tadd, new_k_vals, dt, fcn, ...
+                                   options.funarguments);
+        output_t(iout + (1:iadd)) = tadd;
+        iout = length (output_t);
+      else                         # refine = 1
+        iadd = 1;
+        iout += iadd;
+        output_x(:, iout) = ode_x(:, istep);
+        output_t(iout) = ode_t(istep);
+      end
 
-        ## Call Events function only if a valid result has been found.
-        ## Stop integration if eventbreak is true.
-        if (have_EventFcn)
-          solution.event = ...
-            ode_event_handler (options.Events, t(istep), x(:, istep), [],
-                               options.funarguments{:});
-          if (! isempty (solution.event{1}) && solution.event{1} == 1)
-            t(istep) = solution.event{3}(end);
-            x(:, istep) = solution.event{4}(end, :).';
-            solution.unhandledtermination = false;
-            break;
-          endif
+      ## Call OutputFcn
+      if ((options.haveoutputfunction) && (iadd > 0))
+        xadd = output_x(:, (iout-iadd+1):end);
+        tadd = output_t((iout-iadd+1):end);
+        if (! isempty (options.OutputSel))
+          xadd = xadd(options.OutputSel, :);
         endif
+        stop_solve = feval (options.OutputFcn, tadd, xadd, ...
+                            [], options.funarguments{:});
+        if (stop_solve)
+          solution.unhandledtermination = false;
+          terminal_output = true;
+        endif
+      endif
+      if (terminal_event || terminal_output)
+        break;                     # break from main loop
+      endif
 
-        ## Call OutputFcn only if a valid result has been found.
-        ## Stop integration if function returns true.
-        if (options.haveoutputfunction)
-          cnt = options.Refine + 1;
-          approxtime = linspace (t_old, t_new, cnt);
-          approxvals = interp1 ([t_old, t_new],
-                                [x_old, x_new] .',
-                                approxtime, "linear") .';
-          if (isvector (approxvals) && ! isscalar (approxtime))
-            approxvals = approxvals.';
-          endif
-          if (! isempty (options.OutputSel))
-            approxvals = approxvals(options.OutputSel, :);
-          endif
-          stop_solve = false;
-          for ii = 1:numel (approxtime)
-            stop_solve = feval (options.OutputFcn,
-                                approxtime(ii), approxvals(:, ii), [],
-                                options.funarguments{:});
-            if (stop_solve)
-              break;  # break from inner loop
-            endif
-          endfor
-          if (stop_solve)  # Leave main loop
-            solution.unhandledtermination = false;
-            break;
-          endif
-        endif
-
-      endif
 
       ## move to next time-step
       t_old = t_new;
@@ -303,7 +268,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 (ode_t(end))))
       break;
     endif
 
@@ -313,7 +278,7 @@
   endwhile
 
   ## Check if integration of the ode has been successful
-  if (dir * t(end) < dir * tspan(end))
+  if (dir * ode_t(end) < dir * tspan(end))
     if (solution.unhandledtermination == true)
       warning ("integrate_adaptive:unexpected_termination",
                [" Solving was not successful. ", ...
@@ -322,20 +287,14 @@
                 " This may happen if the stepsize becomes too small. ", ...
                 " Try to reduce the value of 'InitialStep'", ...
                 " and/or 'MaxStep' with the command 'odeset'."],
-               t(end), tspan(end));
-    else
-      warning ("integrate_adaptive: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'."],
-               t(end), tspan(end));
+               ode_t(end), tspan(end));
     endif
   endif
 
   ## Set up return structure
-  solution.t = t(:);
-  solution.x = x.';
+  solution.ode_t = ode_t(:);
+  solution.ode_x = ode_x.';
+  solution.output_t = output_t(:);
+  solution.output_x = output_x.';
 
 endfunction
diff -r dd6b37f67db2 -r 2bd7713ccdb8 scripts/ode/private/ode_event_handler.m
--- a/scripts/ode/private/ode_event_handler.m	Sun Sep 25 06:22:25 2022 -0400
+++ b/scripts/ode/private/ode_event_handler.m	Tue Sep 27 19:48:36 2022 -0400
@@ -115,16 +115,47 @@
         else
           retcell{1} = any (term(idx));     # Stop integration or not
         endif
-        idx = idx(1);  # Use first event found if there are multiple.
-        retcell{2}(evtcnt,1) = idx;
-        ## Calculate the time stamp when the event function returned 0 and
-        ## calculate new values for the integration results, we do both by
-        ## a linear interpolation.
-        tnew = t - evt(idx) * (t - told) / (evt(idx) - evtold(idx));
-        ynew = (y - (t - tnew) * (y - yold) / (t - told)).';
-        retcell{3}(evtcnt,1) = tnew;
-        retcell{4}(evtcnt,:) = ynew;
-        evtcnt += 1;
+        evtcntnew = 1;
+        ## Add all events this step to the output.
+        for idx2 = idx                      # Loop through all values of idx
+          ## Calculate the time stamp when the event function returned 0 and
+          ## calculate new values for the integration results, we do both by
+          ## a linear interpolation.
+          tnew = t - evt(idx2) * (t - told) / (evt(idx2) - evtold(idx2));
+          ynew = (y - (t - tnew) * (y - yold) / (t - told)).';
+          tnews(evtcntnew, 1) = tnew;
+          ynews(evtcntnew, :) = ynew;
+          terms(evtcntnew, 1) = term(idx2);
+          evtcntnew += 1;
+        endfor
+        ## Sort by time of event
+        if length (idx) > 1
+          [tnews, idx_sort] = sort (tnews, "ascend");
+          idxs = idx(idx_sort);
+          ynews = ynews(idx_sort,:);
+          terms = terms(idx_sort);
+        else
+          idxs = idx;
+        endif
+        ## Check for terminal events and remove any events after terminal.
+        ## Any events at same time as first terminal event will be retained.
+        idx3 = find (terms, 1);          # Find first terminal event by time
+        if ! isempty (idx3)
+          t_cutoff = tnews(idx3);
+          ## Last index to return
+          evtcntnew = find (tnews == t_cutoff, 1, "last");
+        else
+          evtcntnew = length (terms);         # Return all indices if no terminal
+        endif
+        idxs = idxs(1:evtcntnew);
+        tnews = tnews(1:evtcntnew);
+
+        ## Populate return values with sorted, clipped values
+        evtcntrange = evtcnt - 1 + (1:evtcntnew);
+        evtcnt += evtcntnew;
+        retcell{2}(evtcntrange, 1) = idxs(:);
+        retcell{3}(evtcntrange, 1) = tnews(:);
+        retcell{4}(evtcntrange, :) = ynews(1:evtcntnew,:);
       endif
 
     endif
