# HG changeset patch
# User Carlo de Falco <carlo.defalco@polimi.it>
# Date 1444588751 -7200
# Node ID 756b052037fb9480866bcb1cb06a9fcf2ef863f8
# Parent  43822bda4f6539e0db0c6e7dee0180a1dc894865
avoid stepping beyond end of thspan in ode solvers

* scripts/ode/private/integrate_adaptive.m: make sure
  not to step beyound end of tspan.

diff -r 43822bda4f65 -r 756b052037fb scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m	Sun Oct 11 19:20:27 2015 +0200
+++ b/scripts/ode/private/integrate_adaptive.m	Sun Oct 11 20:39:11 2015 +0200
@@ -66,17 +66,17 @@
 
   t_new = t_old = t = tspan(1);
   x_new = x_old = x = x0(:);
-  
+
   ## get first initial timestep
   dt = starting_stepsize (order, func, t, x, options.AbsTol,
                           options.RelTol, options.vnormcontrol);
   dt = odeget (options, "InitialStep", dt, "fast_not_empty");
-  
+
   dir = odeget (options, "vdirection", [], "fast");
   dt = dir * min (abs (dt), options.MaxStep);
 
   options.comp = 0.0;
-  
+
   ## Factor multiplying the stepsize guess
   facmin = 0.8;
   facmax = 1.5;
@@ -86,7 +86,7 @@
   if (options.vhaveoutputfunction)
     if (options.vhaveoutputselection)
       solution.vretout = x(options.OutputSel,end);
-    else 
+    else
       solution.vretout = x;
     endif
     feval (options.OutputFcn, tspan, solution.vretout,
@@ -102,14 +102,14 @@
   if (options.vhavenonnegative)
     nn = options.NonNegative;
   endif
-  
+
   solution.vcntloop = 2;
   solution.vcntcycles = 1;
   solution.vcntsave = 2;
   solution.vunhandledtermination = true;
   ireject = 0;
 
-  k_vals = []; 
+  k_vals = [];
   iout = istep = 1;
   while (dir * t_old < dir * tspan(end))
 
@@ -128,20 +128,20 @@
 
     err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol,
                        options.vnormcontrol, x_est);
-    
+
     ## Accepted solution only if err <= 1.0
     if (err <= 1)
 
       solution.vcntloop++;
       ireject = 0;
-            
+
       ## if output time steps are fixed
       if (fixed_times)
 
         t_caught = find ((tspan(iout:end) > t_old)
                          & (tspan(iout:end) <= t_new));
         t_caught = t_caught + iout - 1;
-        
+
         if (! isempty (t_caught))
           t(t_caught) = tspan(t_caught);
           iout = max (t_caught);
@@ -168,7 +168,7 @@
                 t = t(1:id);
                 x(:, id) = solution.vevent{4}(end, :).';
                 x = x(:,1:id);
-                solution.vunhandledtermination = false; 
+                solution.vunhandledtermination = false;
                 break_loop = true;
                 break;
               endif
@@ -198,7 +198,7 @@
           endif
 
         endif
-        
+
       else
 
         t(++istep)  = t_new;
@@ -215,7 +215,7 @@
                   && solution.vevent{1} == 1)
                 t(istep) = solution.vevent{3}(end);
                 x(:, istep) = solution.vevent{4}(end, :).';
-                solution.vunhandledtermination = false; 
+                solution.vunhandledtermination = false;
                 break;
               endif
         endif
@@ -248,7 +248,7 @@
 
       solution.vcntloop = solution.vcntloop + 1;
       vcntiter = 0;
-            
+
     else
 
       ireject++;
@@ -267,14 +267,18 @@
       endif
 
     endif
-    
+
     ## Compute next timestep, formula taken from Hairer
     err += eps;    # avoid divisions by zero
     dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1))));
-    dt = dir * min (abs (dt), options.MaxStep);    
+    dt = dir * min (abs (dt), options.MaxStep);
+
+    ## make sure we don't go past tpan(end)
+    dt = dir * min (abs (dt), abs (tspan(end) - t_old));
 
   endwhile
-  
+
+
   ## Check if integration of the ode has been successful
   if (dir * t(end) < dir * tspan(end))
     if (solution.vunhandledtermination == true)
@@ -300,6 +304,6 @@
   ## Set up return structure
   solution.t = t(:);
   solution.x = x.';
-  
+
 endfunction