# HG changeset patch
# User marek_ka
# Date 1663637242 14400
#      Mon Sep 19 21:27:22 2022 -0400
# Node ID 009e2be83f658a9104e6a41a030fc9efbfc3b317
# Parent  1e25b135b17fec8d4c19a7bdbed0d75f6bccf490
Add Refine to ode tools. Also fixes some Matlab dissimilarities with Events
and OutputFcn results.

diff -r 1e25b135b17f -r 009e2be83f65 scripts/ode/ode23.m
--- a/scripts/ode/ode23.m	Mon Sep 19 17:20:08 2022 -0400
+++ b/scripts/ode/ode23.m	Mon Sep 19 21:27:22 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,11 +260,11 @@
   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
@@ -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
diff -r 1e25b135b17f -r 009e2be83f65 scripts/ode/ode23s.m
--- a/scripts/ode/ode23s.m	Mon Sep 19 17:20:08 2022 -0400
+++ b/scripts/ode/ode23s.m	Mon Sep 19 21:27:22 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,11 +236,11 @@
   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
@@ -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
diff -r 1e25b135b17f -r 009e2be83f65 scripts/ode/ode45.m
--- a/scripts/ode/ode45.m	Mon Sep 19 17:20:08 2022 -0400
+++ b/scripts/ode/ode45.m	Mon Sep 19 21:27:22 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,11 +259,11 @@
   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
@@ -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
diff -r 1e25b135b17f -r 009e2be83f65 scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m	Mon Sep 19 17:20:08 2022 -0400
+++ b/scripts/ode/private/integrate_adaptive.m	Mon Sep 19 21:27:22 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,109 +161,85 @@
 
       solution.cntloop += 1;
       ireject = 0;              # Clear reject counter
-      break_loop = false;
-
-      ## 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 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), [],
+      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;
-              endif
-            endif
-            ## 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 (break_loop)
-              break; # break from inner loop
-            endif
-          endfor
-          
-          if (break_loop)
-            break;   # Leave main loop
-          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_loop = true;
-          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, 1);
         endif
-
-        ## Call OutputFcn only if a valid result has been found.
-        ## Stop integration if function returns true.
-        if (options.haveoutputfunction)
-          xvals = x(:, istep);
-          if (! isempty (options.OutputSel))
-            xvals = xvals(options.OutputSel, 1);
-          endif
-          stop_solve = false;
-          stop_solve = feval (options.OutputFcn,
-                              t(istep), xvals, [],
-                              options.funarguments{:});
-          if (stop_solve)
-            solution.unhandledtermination = false;
-            break_loop = true;
-          endif
+        stop_solve = feval (options.OutputFcn, tadd, xadd, 
+                            [], options.funarguments{:});
+        if (stop_solve)
+          solution.unhandledtermination = false;
+          terminal_output = true;
         endif
-        
-        if (break_loop)  # Leave main loop
-            break;
-        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;
@@ -281,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
 
@@ -291,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. ", ...
@@ -300,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
