## Copyright (C) 2010 Fotios Kasolis <address@hidden>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {} nestfun (@var{f}, @var{init}, @var{n})
## Generate orbit of a discrete time dynamical system.
##
## Let @var{f} be a user-supplied function that defines a dynamical system.
## Given an initial condition @var{init} and a cycle count @var{n}, nestfun
## returns the matrix 
## @code{[init, f(init), f(f(init)), @dots{}, f^n(init)]}, where each element
## is a column vector of length equal to the dimension of the system.
##
## For example:
## 
## @example
## @group
## f = @@(x) 4.0 * x * (1 - x);
## x = nestfun (f, 0.1, 100);
## plot (x, "k", "linewidth", 5);
## @end group
## @end example
## @end deftypefn

function x = nestfun_edited (f, init, n)
  
  if (nargin != 3)
    print_usage ();
  endif

  # Convert function string name to a handle if necessary
  if (ischar (f))
    f = str2func (f, "global");
  endif

  if (! (isscalar (n) && n == fix (n) && n > 0))
    error ("nestfun_edited: N must be an integer greater than 0");
  endif

  x = [init(:), zeros(numel (init), n - 1)];
  
  for count = 1:n-1
    x(:, count + 1) = f (x(:, count))(:);
  endfor

endfunction

%!assert(nestfun_edited(@(x) 1, 1, 10), ones (1,10))
%!assert(nestfun_edited(@(x) 0, -1, 10), [-1, zeros(1,9)])
%!assert(nestfun_edited(@(x) 2*x, 1, 10), 2.^[0:9])

%% Test input validation
%!error nestfun_edited()
%!error nestfun_edited(1)
%!error nestfun_edited(1,2)
%!error nestfun_edited(1,2,3,4)
%!error nestfun_edited(@sin, 1, ones(2,2))
%!error nestfun_edited(@sin, 1, 10.1)
%!error nestfun_edited(@sin, 1, -10)