--- a/scripts/ode/ode45.m	Fri Oct 02 15:07:37 2015 -0400
+++ b/scripts/ode/ode45.m	Sat Oct 03 07:32:50 2015 +0200
@@ -106,8 +106,8 @@
   endif
 
   if (length (vslot) < 2 && ...
-     (isempty (vodeoptions.TimeStepSize) ...
-      || isempty (vodeoptions.TimeStepNumber)))
+      (isempty (vodeoptions.TimeStepSize) ...
+       || isempty (vodeoptions.TimeStepNumber)))
     error ("OdePkg:InvalidArgument", ...
            "second input argument must be a valid vector");
   elseif (vslot(2) == vslot(1))
@@ -251,8 +251,8 @@
   ## option can be set by the user to another value than default value.
   if (isempty (vodeoptions.InitialStep) && strcmp (integrate_func, "adaptive"))
     vodeoptions.InitialStep = vodeoptions.vdirection* ...
-      starting_stepsize (vorder, vfun, vslot(1), vinit, vodeoptions.AbsTol, ...
-                         vodeoptions.RelTol, vodeoptions.vnormcontrol);
+                              starting_stepsize (vorder, vfun, vslot(1), vinit, vodeoptions.AbsTol, ...
+                                                 vodeoptions.RelTol, vodeoptions.vnormcontrol);
     warning ("OdePkg:InvalidArgument", ...
              "option ''InitialStep'' not set, estimated value %f is used", ...
              vodeoptions.InitialStep);
@@ -359,11 +359,11 @@
     if (vmassdependence) ## constant mass matrices have already
       vmass = @(t,x) vodeoptions.Mass (t, x, vodeoptions.vfunarguments{:});
       vfun = @(t,x) vmass (t, x, vodeoptions.vfunarguments{:}) ...
-        \ vfun (t, x, vodeoptions.vfunarguments{:});
+             \ vfun (t, x, vodeoptions.vfunarguments{:});
     else                 ## if (vmassdependence == false)
       vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:});
       vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ...
-        \ vfun (t, x, vodeoptions.vfunarguments{:});
+             \ vfun (t, x, vodeoptions.vfunarguments{:});
     endif
   endif
 
@@ -385,11 +385,11 @@
   ## Postprocessing, do whatever when terminating integration algorithm
   if (vodeoptions.vhaveoutputfunction) ## Cleanup plotter
     feval (vodeoptions.OutputFcn, solution.t(end), ...
-      solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
+           solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
   endif
   if (vodeoptions.vhaveeventfunction)  ## Cleanup event function handling
     odepkg_event_handle (vodeoptions.Events, solution.t(end), ...
-      solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
+                         solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
   endif
 
   ## Print additional information if option Stats is set

--- a/scripts/ode/private/integrate_adaptive.m	Fri Oct 02 15:07:37 2015 -0400
+++ b/scripts/ode/private/integrate_adaptive.m	Sat Oct 03 07:32:50 2015 +0200
@@ -69,7 +69,7 @@
   ## first values for time and solution
   t = tspan(1);
   x = x0(:);
-
+  
   ## get first initial timestep
   dt = odeget (options, "InitialStep",
                starting_stepsize (order, func, t, x, options.AbsTol,
@@ -120,15 +120,20 @@
   z = t;
   u = x;
 
-  k_vals = feval (func, t , x, options.vfunarguments{:});
-
+  k_vals = []; 
+  
   while (counter <= k)
     facmax = 1.5;
 
     ## compute integration step from t to t+dt
-    [s, y, y_est, k_vals] = stepper (func, z(end), u(:,end),
-                                     dt, options, k_vals);
-
+    if (isempty (k_vals))
+      [s, y, y_est, new_k_vals] = stepper (func, z(end), u(:,end),
+                                           dt, options);
+    else
+      [s, y, y_est, new_k_vals] = stepper (func, z(end), u(:,end),
+                                           dt, options, k_vals);
+    endif
+    
     if (options.vhavenonnegative)
       x(options.NonNegative,end) = abs (x(options.NonNegative,end));
       y(options.NonNegative,end) = abs (y(options.NonNegative,end));
@@ -141,27 +146,29 @@
 
     err = AbsRel_Norm (y(:,end), u(:,end), options.AbsTol, options.RelTol,
                        options.vnormcontrol, y_est(:,end));
-
+    
     ## solution accepted only if the error is less or equal to 1.0
     if (err <= 1)
-
+      
       [tk, comp] = kahan (tk, comp, dt);
       options.comp = comp;
       s(end) = tk;
 
       ## values on this interval for time and solution
-      z = [z(end);s];
-      u = [u(:,end),y];
-
+      z = [z(end); s];
+      u = [u(:,end), y];
+      k_vals = new_k_vals;
+      
       ## if next tspan value is caught, update counter
       if ((z(end) == tspan(counter))
           || (abs (z(end) - tspan(counter)) /
               (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) )
         counter++;
-  
-      ## if there is an element in time vector at which the solution is required
-      ## the program must compute this solution before going on with next steps
+        
+        ## if there is an element in time vector at which the solution is required
+        ## the program must compute this solution before going on with next steps
       elseif (vdirection * z(end) > vdirection * tspan(counter))
+
         ## initialize counter for the following cycle
         i = 2;
         while (i <= length (z))
@@ -179,9 +186,9 @@
             ## choose interpolation scheme according to order of the solver
             switch order
               case 1
-               u_interp = linear_interpolation ([z(i-1) z(i)],
-                                                [u(:,i-1) u(:,i)],
-                                                tspan(counter));
+                u_interp = linear_interpolation ([z(i-1) z(i)],
+                                                 [u(:,i-1) u(:,i)],
+                                                 tspan(counter));
               case 2
                 if (! isempty (k_vals))
                   der = k_vals(:,1);
@@ -194,10 +201,10 @@
                                                     der, tspan(counter));
               case 3
                 u_interp = ...
-                  hermite_cubic_interpolation ([z(i-1) z(i)],
-                                               [u(:,i-1) u(:,i)],
-                                               [k_vals(:,1) k_vals(:,end)],
-                                               tspan(counter));
+                hermite_cubic_interpolation ([z(i-1) z(i)],
+                                             [u(:,i-1) u(:,i)],
+                                             [k_vals(:,1) k_vals(:,end)],
+                                             tspan(counter));
               case 4
                 ## if ode45 is used without local extrapolation this function
                 ## doesn't require a new function evaluation.
@@ -217,10 +224,10 @@
                                    - (1776094331/19743644256) * k_vals(:,6)
                                    + (11237099/235043384) * k_vals(:,7));
                 u_interp = ...
-                  hermite_quartic_interpolation ([z(i-1) z(i)],
-                                                 [u(:,i-1) u_half u(:,i)],
-                                                 [k_vals(:,1) k_vals(:,end)],
-                                                 tspan(counter));
+                hermite_quartic_interpolation ([z(i-1) z(i)],
+                                               [u(:,i-1) u_half u(:,i)],
+                                               [k_vals(:,1) k_vals(:,end)],
+                                               tspan(counter));
 
                 ## it is also possible to do a new function evaluation and use
                 ## the quintic hermite interpolator
@@ -239,16 +246,16 @@
                 der(:,2) = feval (func, z(i) , u(:,i),
                                   options.vfunarguments{:});
                 u_interp = ...
-                  hermite_cubic_interpolation ([z(i-1) z(i)],
-                                               [u(:,i-1) u(:,i)],
-                                               der, tspan(counter));
+                hermite_cubic_interpolation ([z(i-1) z(i)],
+                                             [u(:,i-1) u(:,i)],
+                                             der, tspan(counter));
             endswitch
 
             ## add the interpolated value of the solution
             u = [u(:,1:i-1), u_interp, u(:,i:end)];
             
             ## add the time requested
-            z = [z(1:i-1);tspan(counter);z(i:end)];
+            z = [z(1:i-1); tspan(counter); z(i:end)];
 
             ## update counters
             counter++;
@@ -310,7 +317,7 @@
       ## true
       if (options.vhaveeventfunction)
         solution.vevent = odepkg_event_handle (options.Events, t(end),
-            x(:,end), [], options.vfunarguments{:});
+                                               x(:,end), [], options.vfunarguments{:});
         if (! isempty (solution.vevent{1})
             && solution.vevent{1} == 1)
           t(solution.vcntloop-1,:) = solution.vevent{3}(end,:);
@@ -321,13 +328,15 @@
       endif
       
     else
+      
       facmax = 1.0;
+      
     endif
     
     ## Compute next timestep, formula taken from Hairer
     err += eps;    # adding an eps to avoid divisions by zero
-    dt = dt * min (facmax, max (facmin,
-                                fac * (1 / err)^(1 / (order + 1))));
+    dt = dt * min (facmax,
+                   max (facmin, fac * (1 / err)^(1 / (order + 1))));
     dt = vdirection * min (abs (dt), options.MaxStep);
     
     ## Update counters that count the number of iteration cycles
@@ -383,7 +392,3 @@
   solution.x = x(:,1:end-f)';
   
 endfunction
-
-## Local Variables: ***
-## mode: octave ***
-## End: ***

--- a/scripts/ode/private/runge_kutta_45_dorpri.m	Fri Oct 02 15:07:37 2015 -0400
+++ b/scripts/ode/private/runge_kutta_45_dorpri.m	Sat Oct 03 07:32:50 2015 +0200
@@ -63,7 +63,8 @@
 ##
 ## @seealso{odepkg}
 
-function varargout = runge_kutta_45_dorpri (f, t, x, dt, varargin)
+function [t_out, x_out, x_est, k] = ...
+         runge_kutta_45_dorpri (f, t, x, dt, varargin)
 
   persistent a = [0           0          0           0        0          0;
                   1/5         0          0           0        0          0;
@@ -83,17 +84,19 @@
   s = t + dt * b;
   cc = dt * c;
   aa = dt * a;
-  
-  args = varargin{1}.vfunarguments;
   k = zeros (rows (x), 7);
 
-  if (nargin == 5) # only the options are passed
-    k(:,1) = feval (f, t, x, args{:});
-  elseif (nargin == 6) # both the options and the k values are passed
-    k(:,1) = varargin{2}(:,end); # FSAL property
+  if (nargin >= 5) # options are passed
+    args = varargin{1}.vfunarguments;
+    if (nargin >= 6) # both the options and the k values are passed
+      k(:,1) = varargin{2}(:,end); # FSAL property
+    else      
+      k(:,1) = feval (f, t, x, args{:});
+    endif
+  else
+    args = {};
   endif
-  
-  k(:,1) = feval (f, s(1), x, args{:});
+
   k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).', args{:});
   k(:,3) = feval (f, s(3), x + k(:,1:2) * aa(3, 1:2).', args{:});
   k(:,4) = feval (f, s(4), x + k(:,1:3) * aa(4, 1:3).', args{:});
@@ -101,16 +104,15 @@
   k(:,6) = feval (f, s(6), x + k(:,1:5) * aa(6, 1:5).', args{:});
 
   ## compute new time and new values for the unkwnowns
-  varargout{1} = t + dt;
-  varargout{2} = x + k(:,1:6) * cc(:); # 5th order approximation
+  t_out = t + dt;
+  x_out = x + k(:,1:6) * cc(:); # 5th order approximation
 
   ## if the estimation of the error is required
   if (nargout >= 3)
     ## new solution to be compared with the previous one
-    k(:,7) = feval (f, t + dt, varargout{2}, args{:});
+    k(:,7) = feval (f, t + dt, x_out, args{:});
     cc_prime = dt * c_prime;
-    varargout{3} = x + k * cc_prime(:); # x_est
-    varargout{4} = k;
+    x_est = x + k * cc_prime(:); # x_est
   endif
   
 endfunction