#include <octave/oct.h>
#include <octave/parse.h>

DEFUN_DLD(nestfun, args, nargout, "nestfun")
{
  int nargin = args.length();
  
  if (nargin != 3) {
    print_usage ();
  }
  
  octave_idx_type rows = args(1).array_value ().nelem ();
  octave_idx_type cols = args(2).array_value ().elem (0) + 1;
  
  dim_vector dv (2);
  dv(0) = rows; 
  dv(1) = cols;
  
  NDArray x (dv);
  
  dim_vector dvv (2);
  dvv(0) = rows; 
  dvv(1) = 1;
  NDArray v (dvv);
  
  
  for (octave_idx_type i = 0; i <= rows; i++) {
    x(i, 0) = args(1).array_value (). elem(i);
  }
  
  octave_function *fcn = args(0).function_value ();
  octave_value_list fcnargs;
  
  for (octave_idx_type count = 1; count <= cols; count++) {
	
	
    for (octave_idx_type i = 0; i <= rows; i++) {
	  v(i,0) = x(i, count-1);
	}
	
	fcnargs = feval (fcn, octave_value(v), nargout);

	for (octave_idx_type i = 0; i <= rows; i++) {
	  x(i, count) = fcnargs(0).array_value().elem(i,0);
	}
  }
  
  return octave_value(x);
}