FSAL - new stepper implementation

As stated in the previous post, the implementation of the steppers as it was did not allow the possibility to exploit the FSAL (First Same As Last) property of the Bogacki-Shampine algorithm (ode23) and of the Dormand-Prince algorithm (ode45).
The input and output arguments of the steppers have then be modified. As an example here is the runge_kutta_23 stepper:

function varargout = runge_kutta_23 (f, t, x, dt, varargin)
  options = varargin{1};
  k = zeros (size (x, 1), 4);

  if (nargin == 5) % only the options are passed
    k(:,1) = feval (f, t , x, options.vfunarguments{:});
  elseif (nargin == 6) % both the options and the k values are passed
    k(:,1) = varargin{2}(:,end); % FSAL property
  endif
  k(:,2) = feval (f, t + (1/2)*dt, x + dt*(1/2)*k(:,1), options.vfunarguments{:});
  k(:,3) = feval (f, t + (3/4)*dt, x + dt*(3/4)*k(:,2), options.vfunarguments{:});

  %# computing new time and new values for the unkwnowns
  varargout{1} = t + dt; %t_next
  varargout{2} = x + dt.*((2/9)*k(:,1) + (1/3)*k(:,2) + (4/9)*k(:,3)); % return the 3rd order approximation x_next

  %# if the estimation of the error is required
  if (nargout >= 3)
    %# new solution to be compared with the previous one
    k(:,4) = feval (f, t + dt, varargout{2}, options.vfunarguments{:});
    varargout{3} = x + dt.*((7/24)*k(:,1) + (1/4)*k(:,2) + (1/3)*k(:,3) + (1/8)*k(:,4)); %x_est
    varargout{4} = k;
  endif

endfunction

And the call within the solver becomes:

[s, y, y_est, k_vals] = stepper (func, z(end), u(:,end), dt, options, k_vals);

where k_vals has to be initialized for the first iteration as f(t, x).
This implementation will reduce the number of function evaluation for each step.

Furthermore, after some tests in MATLAB, the return values for the solution and the estimate in the  runge_kutta_23 and runge_kutta_45 steppers have been swapped to automatically perform local extrapolation. The MATLAB functions are in fact of order 3 and 5 respectively.