# HG changeset patch
# User marek_ka
# Date 1663622408 14400
#      Mon Sep 19 17:20:08 2022 -0400
# Node ID 1e25b135b17fec8d4c19a7bdbed0d75f6bccf490
# Parent  95fa3182f00ce043808a96a2fcbca584c998f34b
Fix for but #63063. ODE solvers now process multiple Events during one step
and correctly call OutputFcn when there is also an Event (especially a terminal
one). Outputs in struct form are returned in the correct format.

* scripts/ode/private/ode_event_handler.m: Process multiple Events
* scripts/ode/private/integrate_adaptive.m: Process Events and call OutputFcn
     appropriately
* scripts/ode/ode23.m: Struct output formatted properly
* scripts/ode/ode23s.m: Struct output formatted properly
* scripts/ode/ode45.m: Struct output formatted properly

diff -r 95fa3182f00c -r 1e25b135b17f scripts/ode/ode23.m
--- a/scripts/ode/ode23.m	Mon Sep 19 13:14:16 2022 +0200
+++ b/scripts/ode/ode23.m	Mon Sep 19 17:20:08 2022 -0400
@@ -272,9 +272,9 @@
     varargout{1}.y = solution.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 ();
diff -r 95fa3182f00c -r 1e25b135b17f scripts/ode/ode23s.m
--- a/scripts/ode/ode23s.m	Mon Sep 19 13:14:16 2022 +0200
+++ b/scripts/ode/ode23s.m	Mon Sep 19 17:20:08 2022 -0400
@@ -248,9 +248,9 @@
     varargout{1}.y = solution.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 ();
diff -r 95fa3182f00c -r 1e25b135b17f scripts/ode/ode45.m
--- a/scripts/ode/ode45.m	Mon Sep 19 13:14:16 2022 +0200
+++ b/scripts/ode/ode45.m	Mon Sep 19 17:20:08 2022 -0400
@@ -272,9 +272,9 @@
     varargout{1}.y = solution.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 ();
diff -r 95fa3182f00c -r 1e25b135b17f scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m	Mon Sep 19 13:14:16 2022 +0200
+++ b/scripts/ode/private/integrate_adaptive.m	Mon Sep 19 17:20:08 2022 -0400
@@ -150,6 +150,7 @@
 
       solution.cntloop += 1;
       ireject = 0;              # Clear reject counter
+      break_loop = false;
 
       ## if output time steps are fixed
       if (fixed_times)
@@ -165,16 +166,14 @@
             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);
+          ## Call Events or OutputFcn function(s) only if a valid result
+          ## has been found. Stop integration if eventbreak is true.
+          for idenseout = 1:numel (t_caught)
+            id = t_caught(idenseout);
+            td = t(id);
+            ## Check for Events
+            if (have_EventFcn)
               solution.event = ...
                 ode_event_handler (options.Events, t(id), x(:, id), [],
                                    options.funarguments{:});
@@ -185,41 +184,28 @@
                 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.';
+            ## Call OutputFcn
+            if (options.haveoutputfunction)
+              xvals = x(:,id);
+              if (! isempty (options.OutputSel))
+                xvals = xvals(options.OutputSel, 1);
+              endif
+              stop_solve = feval (options.OutputFcn, t(id), xvals, 
+                                  [], options.funarguments{:});
+              if (stop_solve)
+                solution.unhandledtermination = false;
+                break_loop = true;
+              endif
             endif
-            if (! isempty (options.OutputSel))
-              approxvals = approxvals(options.OutputSel, :);
+            if (break_loop)
+              break; # break from inner loop
             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
+          endfor
+          
+          if (break_loop)
+            break;   # Leave main loop
           endif
 
         endif
@@ -240,37 +226,29 @@
             t(istep) = solution.event{3}(end);
             x(:, istep) = solution.event{4}(end, :).';
             solution.unhandledtermination = false;
-            break;
+            break_loop = 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_new],
-                                [x_old, x_new] .',
-                                approxtime, "linear") .';
-          if (isvector (approxvals) && ! isscalar (approxtime))
-            approxvals = approxvals.';
-          endif
+          xvals = x(:, istep);
           if (! isempty (options.OutputSel))
-            approxvals = approxvals(options.OutputSel, :);
+            xvals = xvals(options.OutputSel, 1);
           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
+          stop_solve = feval (options.OutputFcn,
+                              t(istep), xvals, [],
+                              options.funarguments{:});
+          if (stop_solve)
             solution.unhandledtermination = false;
+            break_loop = true;
+          endif
+        endif
+        
+        if (break_loop)  # Leave main loop
             break;
-          endif
         endif
 
       endif
diff -r 95fa3182f00c -r 1e25b135b17f scripts/ode/private/ode_event_handler.m
--- a/scripts/ode/private/ode_event_handler.m	Mon Sep 19 13:14:16 2022 +0200
+++ b/scripts/ode/private/ode_event_handler.m	Mon Sep 19 17:20:08 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
