# HG changeset patch
# User marek_ka
# Date 1663801203 14400
#      Wed Sep 21 19:00:03 2022 -0400
# Node ID ee90d0e70d5042e8e41ecc497877fff8fa1a9084
# Parent  a103cfd738c5e7e78df614b162028e3bc208e04c
Bug fix for #63063 (including #49408)

diff -r a103cfd738c5 -r ee90d0e70d50 scripts/ode/ode23.m
--- a/scripts/ode/ode23.m	Wed Sep 21 13:48:27 2022 -0700
+++ b/scripts/ode/ode23.m	Wed Sep 21 19:00:03 2022 -0400
@@ -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,8 @@
     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}.x = solution.ode_t.';  # Time stamps saved in field x (row vector)
+    varargout{1}.y = solution.ode_x.';  # Results are saved in field y (row vector)
     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"))
@@ -370,6 +366,12 @@
 %!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
+%! assert ( (length (t1) == length (t3)) && ...
+%!          (length (t2) == 4*length (t1) - 3) );
 %!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 +437,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,10 +454,9 @@
 %! 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");
diff -r a103cfd738c5 -r ee90d0e70d50 scripts/ode/ode23s.m
--- a/scripts/ode/ode23s.m	Wed Sep 21 13:48:27 2022 -0700
+++ b/scripts/ode/ode23s.m	Wed Sep 21 19:00:03 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,8 @@
     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 +236,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 vector)
+    varargout{1}.y = solution.ode_x.';  # Results are saved in field y (row vector)
     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 +258,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 +384,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 +400,12 @@
 %!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
+%! assert ( (length (t1) == length (t3)) && ...
+%!          (length (t2) == 4*length (t1) - 3) );
 %!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 +437,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");
diff -r a103cfd738c5 -r ee90d0e70d50 scripts/ode/ode45.m
--- a/scripts/ode/ode45.m	Wed Sep 21 13:48:27 2022 -0700
+++ b/scripts/ode/ode45.m	Wed Sep 21 19:00:03 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 = ...
@@ -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,8 @@
     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 +259,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 vector)
+    varargout{1}.y = solution.ode_x.';  # Results are saved in field y (row vector)
     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 +281,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 +349,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 +366,14 @@
 %! [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
+%! assert ( (length (t2) == length (t3)) && ...
+%!          (length (t2) == 4*length (t1) - 3) );
 %!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 +393,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 +448,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,10 +465,9 @@
 %! 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");
@@ -511,7 +511,7 @@
 %!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)));
 
 %!error <Invalid call> ode45 ()
diff -r a103cfd738c5 -r ee90d0e70d50 scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m	Wed Sep 21 13:48:27 2022 -0700
+++ b/scripts/ode/private/integrate_adaptive.m	Wed Sep 21 19:00:03 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);
@@ -94,13 +94,24 @@
   facmin = 0.8;
   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,131 +161,85 @@
 
       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;
-
-        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;
-
-          ## 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), [],
+      terminal_event = false;
+      terminal_output = false;
+      
+      istep++;
+      ode_t(istep) = t_new;
+      ode_x(:, istep) = x_new;
+      iadd = 0;                 # Number of output points added this iteration
+      
+      ## Check for Events
+      if (have_EventFcn)
+        solution.event = ...
+          ode_event_handler (options.Events, t_new, x_new, [],
                                    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
-
-          ## 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
-
+        ## 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
+      
+      ## 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
+        ## 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
-
-      else   # not fixed times
-
-        t(++istep)  = t_new;
-        x(:, istep) = x_new;
-        iout = istep;
-
-        ## 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
+      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 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
-
-        ## 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
+        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
 
+        
       ## move to next time-step
       t_old = t_new;
       x_old = x_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 a103cfd738c5 -r ee90d0e70d50 scripts/ode/private/ode_event_handler.m
--- a/scripts/ode/private/ode_event_handler.m	Wed Sep 21 13:48:27 2022 -0700
+++ b/scripts/ode/private/ode_event_handler.m	Wed Sep 21 19:00:03 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
