function [ t, y, schemeData ] = ...
       odeCFLmultipleSteps(intFunc, schemeFunc, tspan, y0, options, schemeData)
% odeCFLmultipleSteps: Handle the length(tspan) > 2 case for odeCFLn.
%
% [ t, y, schemeData ] = ...
%      odeCFLmultipleSteps(intFunc, schemeFunc, tspan, y0, options, schemeData)
%
% Makes repeated calls to an odeCFLn routine to integrate a system forward
%   in time, stopping at specified times to record the data.
%
% This routine is not intended to be called directly by the toolbox user;
%   rather, it is some common code for the odeCFLn routines that was factored
%   into a separate routine.
%
% Apart from the parameter intFunc, the parameters for this routine are
%   the same as those for odeCFLn.
%
% parameters:
%   intFunc      Function handle to an odeCFLn integration routine for
%                  CFL constrained ODE systems.  This routine is used
%                  to integrate between the pairs of times in tspan.
%   schemeFunc	 Function handle to a CFL constrained ODE system
%                  (typically an approximation to an HJ term, see below).
%   tspan        Range of time over which to integrate (see below).
%   y0           Initial condition vector 
%                  (typically the data array in vector form).
%   options      An option structure generated by odeCFLset 
%                  (use [] as a placeholder if necessary).
%   schemeData   Structure passed through to schemeFunc.
%
%
%   t            Output time(s) (see below).
%   y            Output state (see below).
%   schemeData   Output version of schemeData (see below).
%
% The time interval tspan must be given as
%   A vector with three or more entries, in which case the output will
%   be column vector t = tspan and each row of y will be the solution
%   at one of the times in tspan.  Unlike Matlab's ode suite routines,
%   odeCFLmultipleSteps just repeatedly calls the function handle intfunc
%   to accomplish this goal, so it is not particularly efficient.

% Copyright 2005 Ian M. Mitchell (mitchell@cs.ubc.ca).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Factored from odeCFLn, Ian Mitchell, 12/06/04.

  %---------------------------------------------------------------------------
  % Number of timesteps to be returned.
  numT = length(tspan);

  %---------------------------------------------------------------------------
  % If we were asked for the solution at multiple timesteps,
  %   call back for each pair of timesteps.
  if(numT > 2)

    t = reshape(tspan, numT, 1);

    if(iscell(y))
      numY = length(y);
      y = cell(numY, 1);
      for i = 1 : numY
        y{i} = zeros(numT, length(y0{i}));
        y{i}(1,:) = y0{i}';
      end
    else
      y = zeros(numT, length(y0));
      y(1,:) = y0';
    end

    yout = y;
    for n = 2 : numT
      [ t(n), yout, schemeData ] = ...
       feval(intFunc, schemeFunc, [ t(n-1), t(n) ], yout, schemeData, options);

      if(iscell(y))
        for i = 1 : numY
          y{i}(n,:) = yout{i}';
        end
      else
        y(n,:) = yout';
      end

    end
  
  %---------------------------------------------------------------------------
  else
    % This routine is only for finding the solution at multiple timesteps.
    error('tspan must contain at least three entries');
  end
