# HG changeset patch
# User Ken Marek <marek_ka@mercer.edu>
# Date 1664997827 14400
#      Wed Oct 05 15:23:47 2022 -0400
# Node ID ea8d7a829121e08c5e6f30e1a32cded8ac048525
# Parent  f4b71c459b88597a49f3bd47d7d765d230fec7d3
Preliminary bug patch, bug #63162

diff -r f4b71c459b88 -r ea8d7a829121 scripts/ode/ode23.m
--- a/scripts/ode/ode23.m	Tue Sep 27 19:48:36 2022 -0400
+++ b/scripts/ode/ode23.m	Wed Oct 05 15:23:47 2022 -0400
@@ -101,6 +101,11 @@
   solver = "ode23";
   order  = 3;
 
+  ## FIXME: Adding additional function arguments is not documented, or
+  ## compatible with Matlab documentation (I don't know about implementation).
+  ## Using additional arguments should also be possible by constructing
+  ## an anonymous function, which should probably be the preferred
+  ## method.
   if (nargin >= 4)
     if (! isstruct (varargin{1}))
       ## varargin{1:len} are parameters for fcn
@@ -238,12 +243,10 @@
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
+    feval (odeopts.OutputFcn, [], [], "done");
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
-    ode_event_handler (odeopts.Events, solution.ode_t(end), ...
-                       solution.ode_x(end,:).', "done", ...
-                       odeopts.funarguments{:});
+    ode_event_handler ([], [], [], [], [], "done");
   endif
 
   ## Print additional information if option Stats is set
@@ -460,11 +463,11 @@
 %! opt = odeset ("Events", @fevn, "NormControl", "on");
 %! sol = ode23 (@fpol, [0 10], [2 0], opt);
 %! assert ([sol.ie, sol.xe, sol.ye.'],
-%!         [2.0, 2.496110, -0.830550, -2.677589], .5e-1);
+%!         [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
 %!test  # Events option, five output arguments
 %! 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);
+%! assert ([vie, vxe, ye], [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
 %!test  # Mass option as function
 %! opt = odeset ("Mass", @fmas);
 %! sol = ode23 (@fpol, [0 2], [2 0], opt);
diff -r f4b71c459b88 -r ea8d7a829121 scripts/ode/ode23s.m
--- a/scripts/ode/ode23s.m	Tue Sep 27 19:48:36 2022 -0400
+++ b/scripts/ode/ode23s.m	Wed Oct 05 15:23:47 2022 -0400
@@ -102,6 +102,11 @@
   solver = "ode23s";
   order = 2;
 
+  ## FIXME: Adding additional function arguments is not documented, or
+  ## compatible with Matlab documentation (I don't know about implementation).
+  ## Using additional arguments should also be possible by constructing
+  ## an anonymous function, which should probably be the preferred
+  ## method.
   if (nargin >= 4)
     if (! isstruct (varargin{1}))
       ## varargin{1:len} are parameters for fcn
@@ -214,12 +219,10 @@
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
+    feval (odeopts.OutputFcn, [], [], "done");
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
-    ode_event_handler (odeopts.Events, solution.ode_t(end), ...
-                       solution.ode_x(end,:).', "done", ...
-                       odeopts.funarguments{:});
+    ode_event_handler ([], [], [], [], [], "done");
   endif
 
   ## Print additional information if option Stats is set
@@ -456,6 +459,15 @@
 %! assert (sol.ie(1), 2);
 %! assert (isfield (sol, "xe"));
 %! assert (isfield (sol, "ye"));
+%!test  # Events option, now stop integration
+%! opt = odeset ("Events", @fevn, "NormControl", "on");
+%! sol = ode23s (@fpol, [0 10], [2 0], opt);
+%! assert ([sol.ie, sol.xe, sol.ye.'],
+%!         [2.0, 9.094439, -0.996480, -14.180147], 2e-2);
+%!test  # Events option, five output arguments
+%! opt = odeset ("Events", @fevn, "NormControl", "on");
+%! [t, y, vxe, ye, vie] = ode23s (@fpol, [0 10], [2 0], opt);
+%! assert ([vie, vxe, ye], [2.0, 9.094439, -0.996480, -14.180147], 2e-2);
 %!test  # Mass option as function
 %! opt = odeset ("Mass", @fmas);
 %! sol = ode23s (@fpol, [0 2], [2 0], opt);
diff -r f4b71c459b88 -r ea8d7a829121 scripts/ode/ode45.m
--- a/scripts/ode/ode45.m	Tue Sep 27 19:48:36 2022 -0400
+++ b/scripts/ode/ode45.m	Wed Oct 05 15:23:47 2022 -0400
@@ -98,6 +98,11 @@
   solver = "ode45";
   order  = 5;  # runge_kutta_45_dorpri uses local extrapolation
 
+  ## FIXME: Adding additional function arguments is not documented, or
+  ## compatible with Matlab documentation (I don't know about implementation).
+  ## Using additional arguments should also be possible by constructing
+  ## an anonymous function, which should probably be the preferred
+  ## method.
   if (nargin >= 4)
     if (! isstruct (varargin{1}))
       ## varargin{1:len} are parameters for fcn
@@ -237,12 +242,10 @@
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
+    feval (odeopts.OutputFcn, [], [], "done");
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
-    ode_event_handler (odeopts.Events, solution.ode_t(end), ...
-                       solution.ode_x(end,:).', "done", ...
-                       odeopts.funarguments{:});
+    ode_event_handler ([], [], [], [], [], "done");
   endif
 
   ## Print additional information if option Stats is set
@@ -471,12 +474,12 @@
 %! opt = odeset ("Events", @fevn, "NormControl", "on");
 %! sol = ode45 (@fpol, [0 10], [2 0], opt);
 %! assert ([sol.ie, sol.xe, sol.ye.'],
-%!         [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
+%!         [2.0, 2.496110, -0.830550, -2.677589], 2e-3);
 %!test  # Events option, five output arguments
 %! opt = odeset ("Events", @fevn, "NormControl", "on");
 %! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt);
 %! assert ([vie, vxe, ye],
-%!         [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
+%!         [2.0, 2.496110, -0.830550, -2.677589], 2e-3);
 %!test  # Mass option as function
 %! opt = odeset ("Mass", @fmas);
 %! sol = ode45 (@fpol, [0 2], [2 0], opt);
diff -r f4b71c459b88 -r ea8d7a829121 scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m	Tue Sep 27 19:48:36 2022 -0400
+++ b/scripts/ode/private/integrate_adaptive.m	Wed Oct 05 15:23:47 2022 -0400
@@ -113,16 +113,14 @@
     else
       solution.retout = output_x;
     endif
-    feval (options.OutputFcn, tspan, solution.retout, "init",
-           options.funarguments{:});
+    feval (options.OutputFcn, tspan, solution.retout, "init");
   endif
 
   ## Initialize the EventFcn
   have_EventFcn = false;
   if (! isempty (options.Events))
     have_EventFcn  = true;
-    ode_event_handler (options.Events, tspan(1), ode_x,
-                       "init", options.funarguments{:});
+    ode_event_handler (tspan(1), ode_x, [], order, options.Events, "init");
   endif
 
   if (options.havenonnegative)
@@ -171,8 +169,8 @@
 
       ## Check for Events
       if (have_EventFcn)
-        solution.event = ode_event_handler (options.Events, t_new, x_new, ...
-                                            [], options.funarguments{:});
+        solution.event = ode_event_handler (t_new, x_new, new_k_vals, ...
+                                            [], []);
         ## Check for terminal Event
         if (! isempty (solution.event{1}) && solution.event{1} == 1)
           ode_t(istep) = solution.event{3}(end);
@@ -193,8 +191,7 @@
           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);
+                                     output_t(t_caught), new_k_vals);
         endif
         ## Add a possible additional output value if we found a terminal Event
         if ((terminal_event == true) && ...
@@ -210,8 +207,7 @@
         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);
+                                   tadd, new_k_vals);
         output_t(iout + (1:iadd)) = tadd;
         iout = length (output_t);
       else                         # refine = 1
@@ -228,8 +224,7 @@
         if (! isempty (options.OutputSel))
           xadd = xadd(options.OutputSel, :);
         endif
-        stop_solve = feval (options.OutputFcn, tadd, xadd, ...
-                            [], options.funarguments{:});
+        stop_solve = feval (options.OutputFcn, tadd, xadd, []);
         if (stop_solve)
           solution.unhandledtermination = false;
           terminal_output = true;
diff -r f4b71c459b88 -r ea8d7a829121 scripts/ode/private/ode_event_handler.m
--- a/scripts/ode/private/ode_event_handler.m	Tue Sep 27 19:48:36 2022 -0400
+++ b/scripts/ode/private/ode_event_handler.m	Wed Oct 05 15:23:47 2022 -0400
@@ -24,43 +24,45 @@
 ########################################################################
 
 ## -*- texinfo -*-
-## @deftypefn {} {@var{retval} =} ode_event_handler (@var{@@evtfcn}, @var{t}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{})
+## @deftypefn {} {@var{retval} =} ode_event_handler (@var{t}, @var{y}, @var{k_vals}, @var{ord}, @var{@@efcn}, @var{flag})
 ##
-## Return the solution of the event function (@var{@@evtfcn}) which is
+## Return the solution of the event function (@var{@@efcn}) which is 
 ## specified in the form of a function handle.
 ##
-## The second input argument @var{t} is a scalar double and specifies the time
+## The first input argument @var{t} is a scalar double and specifies the time
 ## of the event evaluation.
 ##
-## 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 second input argument @var{y} is a scalar or a column vector of type
+## double which specifies the solution(s) at time @var{t}.
 ##
-## The fourth input argument @var{flag} is of type string.  Valid values are:
+## The third input argument @var{k_vals} is a vector or matrix with the
+## k values obtained from the most recent integration step.
+##
+## The fourth input argument @var{ord} the order of the integration technique.
+##
+## The sixth input argument @var{flag} is of type string.  Valid values are:
 ##
 ## @table @option
 ## @item  @qcode{"init"}
 ## Initialize internal persistent variables of the function
 ## @code{ode_event_handler} and return an empty cell array of size 4.
 ##
-## @item  @qcode{"calc"}
-## Evaluate the event function and return the solution @var{retval} as a cell
-## array of size 4.
+## @item  @qcode{""}
+## (default) Evaluate the event function and return the solution @var{retval}
+## as a cell array of size 4. Inputs @var{ord} and @var{@@efcn} are ignored,
+## since they are set in the @qcode{"init"} step.
 ##
 ## @item  @qcode{"done"}
 ## Clean up internal variables of the function @code{ode_event_handler} and
-## return an empty cell array of size 4.
+## return an empty cell array of size 4. All other inputs are ignored with
+## this flag.
 ## @end table
 ##
-## 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 and it should never be
 ## necessary to call it directly.
 ## @end deftypefn
 
-function retval = ode_event_handler (evtfcn, t, y, flag = "", varargin)
+function retval = ode_event_handler (t, y, k_vals, ord, efcn, flag = "")
 
   ## No error handling has been implemented in this function to achieve
   ## the highest performance possible.
@@ -77,23 +79,14 @@
   ## 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;
+  persistent evtold told yold retcell order evtfcn;
   persistent evtcnt = 1;   # Don't remove.  Required for Octave parser.
   persistent firstrun = true;
 
   if (isempty (flag))
     ## Process the event, i.e.,
     ## find the zero crossings for either a rising or falling edge
-    if (! iscell (y))
-      inpargs = {evtfcn, t, y};
-    else
-      inpargs = {evtfcn, t, y{1}, y{2}};
-      y = y{1};  # Delete cell element 2
-    endif
-    if (nargin > 4)
-      inpargs = {inpargs{:}, varargin{:}};
-    endif
-    [evt, term, dir] = feval (inpargs{:});
+    [evt, term, dir] = evtfcn (t, y);
 
     ## We require that all return values be row vectors
     evt = evt(:).'; term = term(:).'; dir = dir(:).';
@@ -119,10 +112,21 @@
         ## 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)).';
+          ## calculate new values for the integration results. We do both by
+          ## root solving, calling the Event function with y values from
+          ## the RK interpolation polynomial. Set tolerance to zero (actually
+          ## uses machine tolerance based criterion in this case) since we
+          ## don't have a way for the user to specify an acceptable tolerance.
+          ## For simple Event functions, this means we're basically root
+          ## finding on a small interval for a polynomial, so we expect
+          ## pretty quick convergence.
+          tvals = [told t];
+          yvals = [yold y];
+          tnew = fzero(@(t2) evtfcn_val (evtfcn, t2, ...
+                       runge_kutta_interpolate (order, tvals, yvals, ...
+                       t2, k_vals), idx2), tvals, optimset ("TolX", 0));
+          ynew = runge_kutta_interpolate (order, tvals, yvals, tnew, k_vals);
+          
           tnews(evtcntnew, 1) = tnew;
           ynews(evtcntnew, :) = ynew;
           terms(evtcntnew, 1) = term(idx2);
@@ -170,17 +174,10 @@
     ## a value for evtold.
 
     firstrun = true;
+    order = ord;
+    evtfcn = efcn;
 
-    if (! iscell (y))
-      inpargs = {evtfcn, t, y};
-    else
-      inpargs = {evtfcn, t, y{1}, y{2}};
-      y = y{1};  # Delete cell element 2
-    endif
-    if (nargin > 4)
-      inpargs = {inpargs{:}, varargin{:}};
-    endif
-    [evtold, ~, ~] = feval (inpargs{:});
+    [evtold, ~, ~] = evtfcn (t, y);
 
     ## We require that all return values be row vectors
     evtold = evtold(:).'; told = t; yold = y;
@@ -190,9 +187,15 @@
   elseif (strcmp (flag, "done"))
     ## Clear this event handling function
     firstrun = true;
-    evtold = told = yold = evtcnt = [];
+    evtold = told = yold = evtcnt = order = evtfcn = [];
     retval = retcell = cell (1,4);
 
   endif
 
 endfunction
+
+function val = evtfcn_val (evtfcn, t, y, ind)
+  [evt, ~, ~] = evtfcn (t, y);
+  val = evt(ind);
+endfunction
+  
diff -r f4b71c459b88 -r ea8d7a829121 scripts/ode/private/runge_kutta_interpolate.m
--- a/scripts/ode/private/runge_kutta_interpolate.m	Tue Sep 27 19:48:36 2022 -0400
+++ b/scripts/ode/private/runge_kutta_interpolate.m	Wed Oct 05 15:23:47 2022 -0400
@@ -23,22 +23,21 @@
 ##
 ########################################################################
 
-function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, fcn, args)
+function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals)
 
   switch (order)
 
     case 1
+      ## Unused
       u_interp = interp1 (z, u.', t, "linear");
 
     case 2
-      if (! isempty (k_vals))
-        der = k_vals(:,1);
-      else
-        der = feval (fcn, z(1) , u(:,1), args);
-      endif
+      ## ode23s with Rosenbrock scheme:
+      der = k_vals(:,1);
       u_interp = quadratic_interpolation (z, u, der, t);
 
     case 3
+      ## ode23 with Bogacki-Shampine scheme:
       u_interp = hermite_cubic_interpolation (z, u, k_vals, t);
 
     case 5
@@ -46,11 +45,10 @@
       u_interp = hermite_quartic_interpolation (z, u, k_vals, t);
 
     otherwise
-      warning (["High order interpolation not yet implemented: ", ...
-                "using cubic interpolation instead"]);
-      der(:,1) = feval (fcn, z(1), u(:,1), args);
-      der(:,2) = feval (fcn, z(2), u(:,2), args);
-      u_interp = hermite_cubic_interpolation (z, u, der, t);
+      ## This should never happen
+      warning (["Invalid/unimplemented interpolation order: ", ...
+                "using linear interpolation instead"]);
+      u_interp = interp1 (z, u.', t, "linear");
 
   endswitch
 
