diff --git a/pararrayfun.m b/pararrayfun.m
--- a/pararrayfun.m
+++ b/pararrayfun.m
@@ -39,31 +39,91 @@
   endif
 
   [nargs, uniform_output, error_handler, ...
-  verbose_level, chunks_per_proc, vectorized] = parcellfun_opts (varargin); 
+    verbose_level, chunks_per_proc, vectorized, ...
+    cat_dimensions, idx_dimensions] = parcellfun_opts (varargin); 
 
   args = varargin(1:nargs);
   opts = varargin(nargs+1:end);
+  
+  # use default indexing ?
+  default_indexing = isempty(idx_dimensions);
+  
   if (nargs == 0)
     print_usage ();
-  elseif (nargs > 1)
-    [err, args{:}] = common_size (args{:});
-    if (err)
-      error ("pararrayfun: arguments size must match");
+  elseif default_indexing
+    if (nargs > 1)
+      [err, args{:}] = common_size (args{:});
+      if (err)
+        error ("pararrayfun: arguments size must match or\n\
+                    need IdxDimensions option");
+      endif
     endif
+    indexed_args_mask = ones(size(args), "logical");
+    njobs = numel (args{1});
+  else
+    assert(numel(idx_dimensions), numel(args));
+    # indexes of the indexed arguments
+    indexed_args_mask = (idx_dimensions != 0);
+    # sizes of the indexed arguments along the indexed dimension
+    indexed_dimensions_sizes = cellfun(@(a, dim) size(a, dim), 
+                                       args(indexed_args_mask), 
+                                       num2cell(idx_dimensions(indexed_args_mask)),
+                                       "UniformOutput", false);
+    # all these lengths must be equal
+    if numel(indexed_dimensions_sizes) > 1
+      assert(isequal(indexed_dimensions_sizes{:}))
+    endif
+    # number of "jobs" = length along the indexed dimension
+    # they are all equal, just take the first one
+    njobs = indexed_dimensions_sizes{1};
   endif
 
-  njobs = numel (args{1});
+  
 
   if (vectorized && chunks_per_proc > 0 && chunks_per_proc < njobs / nproc)
     ## If "Vectorized" is on, we apply the function directly on chunks of
     ## arrays.
     [varargout{1:nargout}] = chunk_parcellfun (nproc, chunks_per_proc, ...
-    func, error_handler, verbose_level, args{:});
+    func, error_handler, verbose_level, cat_dimensions, idx_dimensions, args{:});
   else
-    args = cellfun (@num2cell, args, "UniformOutput", false,
-    "ErrorHandler",  @arg_class_error);
+    if default_indexing
+      args = cellfun (@num2cell, args, "UniformOutput", false,
+                      "ErrorHandler",  @arg_class_error);
+    else
+      # ## set the split sizes
+      # # FIXME: duplicated code with chunk_parcellfun
+      # split_sizes = cellfun(@(argument) num2cell(size(argument)), 
+                             # args, "UniformOutput", false);
+      # for arg_idx = 1:numel(args)
+        # if (idx_dimensions(arg_idx) != 0)
+          # split_sizes{arg_idx}{idx_dimensions(arg_idx)} = ones(1, njobs);
+        # endif
+      # endfor
+      # # split only indexed arguments (in one element chunks along the indexed dimension)
+      # args(indexed_args_mask) = cellfun (
+             # @(argument, split_size) mat2cell (argument, split_size{:}), ...
+               # args(indexed_args_mask), split_sizes(indexed_args_mask), ...
+               # "UniformOutput", false)
+      # apply num2cell only on indexed arguments
+      args(indexed_args_mask) = cellfun (@num2cell,
+                                 args(indexed_args_mask), "UniformOutput", false);
+    endif
+    
+    # # need the same opts, but without IdxDimensions
+    # opts_parcellfun = {};
+    opts_parcellfun = opts;
+    # for idx_opt = 1:2:numel(opts);
+      # opt = opts{idx_opt};
+      # val = opts{idx_opt+1};
+      # if (strcmp(tolower(opt), "idxdimensions"))
+        # opts_parcellfun = [opts_parcellfun, "IdxDimensions", ones(1, nargs)];
+      # else
+        # opts_parcellfun = [opts_parcellfun, opt, val];
+      # endif
+    # endfor
 
-    [varargout{1:nargout}] = parcellfun (nproc, func, args{:}, opts{:});
+    # call parcellfun
+    [varargout{1:nargout}] = parcellfun (nproc, func, args{:}, opts_parcellfun{:});
   endif
 
 endfunction
@@ -71,3 +131,139 @@
 function arg_class_error (S, X)
   error ("arrayfun: invalid argument of class %s", class (X))
 endfunction
+
+%!test % standard call
+%! f = @(x, y) x + y;
+%! v = 1:5;
+%! r = pararrayfun(2, f, v, v);
+%! assert(r, 2:2:10)
+
+%!function [x, y] = f(z)
+%!  x = z;
+%!  y = z.^2;
+%!endfunction
+%!
+%!test % standard 2 outputs
+%! v = 1:10;
+%! [x, y] = pararrayfun(2, @f, v);
+%! assert(x, v)
+%! assert(y, v.^2)
+%!
+%!test % vectorized, 2 outputs
+%! v = 1:10;
+%! [x, y] = pararrayfun(2, @f, v, "Vectorized", true, "ChunksPerProc", 2);
+%! assert(x, v)
+%! assert(y, v.^2)
+
+%!function [x, y] = f(z)
+%!   x = repmat(z(:), [1, 2, 2]); y = repmat(z(:)', [2, 1, 2]);
+%!endfunction
+%!
+%!test % setting concatenation dimensions for output arguments
+%! [x, y] = pararrayfun(2, @f, [1:6], "Vectorized", true, "ChunksPerProc", 2, "VerboseLevel", 1, "CatDimensions", [1 2]);
+
+%!test % IdxDimensions - not Vectorized
+%! f = @(x, y) sum(x + y.', 2);
+%! m = [1:5; 2:6];
+%! r = pararrayfun(2, f, m.', m, "IdxDimensions", [1 2]);
+%! assert(r, f(m.', m))
+
+%!test % IdxDimensions - Vectorized
+%! f = @(x, y) sum(x + y.', 2);
+%! m = [1:10; 2:11];
+%! r = pararrayfun(2, f, m.', m, "IdxDimensions", [1 2], "ChunksPerProc", 2, "Vectorized", true);
+%! assert(r, (6:4:42).')
+
+%!test % unindexed arguments - Vectorized
+%! f = @(x, y) sum(x + y.', 2);
+%! m = [1:10; 2:11];
+%! s = 2;
+%! r = pararrayfun(2, f, m.', s, "IdxDimensions", [1 0], "ChunksPerProc", 2, "Vectorized", true);
+%! assert(r, (7:2:25).')
+
+%!function [x, y] = f(u, v)
+%!   x = u + v.';
+%!   y = x.';
+%!endfunction
+%!
+%!test % setting CatDimensions and IdxDimensions
+%! u = [1:10; 2:11];
+%! v = u.';
+%! [x, y] = pararrayfun(2, @f, u,v, "VerboseLevel", 1, "CatDimensions", [2 1],"IdxDimensions", [2 1] );
+%! assert(x, [2:2:20; 4:2:22])
+%! assert(y, [2:2:20; 4:2:22].')
+
+%!function [x, y] = f(u, v)
+%!   x = u + v.';
+%!   y = x.';
+%!endfunction
+%!
+%!test % setting CatDimensions and IdxDimensions - vectorized
+%! u = [1:10; 2:11];
+%! v = u.';
+%! [x, y] = pararrayfun(2, @f, u,v, "Vectorized", true, "ChunksPerProc", 2, "VerboseLevel", 1, "CatDimensions", [2 1],"IdxDimensions", [2 1] );
+%! assert(x, [2:2:20; 4:2:22])
+%! assert(y, [2:2:20; 4:2:22].')
+
+%!function [x, y] = f(u, v)
+%!   if isargout(1)
+%!     x = u;
+%!   endif
+%!   if isargout(2)
+%!     y = v;
+%!   endif
+%!endfunction
+%!
+%!test % ignored outputs
+%! u = (1:10).';
+%! v = u;
+%! [~, y] = pararrayfun(2, @f, u,v, "Vectorized", true, "ChunksPerProc", 2, "VerboseLevel", 1);
+%! assert(y, v)
+
+%!function [x, y] = f(u, v)
+%!   if isargout(1)
+%!     x = u + v.';
+%!   endif
+%!   if isargout(2)
+%!     y = x.';
+%!   endif
+%!endfunction
+%!
+%!test % ignored outputs - setting CatDimensions and IdxDimensions
+%! u = [1:10; 2:11];
+%! v = u.';
+%! [~, y] = pararrayfun(2, @f, u,v, "Vectorized", true, "ChunksPerProc", 2, "VerboseLevel", 1, "CatDimensions", [2 1],"IdxDimensions", [2 1] );
+# %! assert(x, [2:2:20; 4:2:22])
+%! assert(y, [2:2:20; 4:2:22].')
+
+%!function [x, y] = f(u, v, w)
+%!   x = u + v.' + w;
+%!   y = x.';
+%!endfunction
+%!
+%!test % setting CatDimensions and IdxDimensions, with non-indexed arg
+%! u = [1:10; 2:11];
+%! v = u.';
+%! w = 3;
+%! [x, y] = pararrayfun(2, @f, u, v, w, "Vectorized", true, "ChunksPerProc", 2, "VerboseLevel", 1, "CatDimensions", [2 1],"IdxDimensions", [2 1 0] );
+%! assert(x, [2:2:20; 4:2:22]+3)
+%! assert(y, [2:2:20; 4:2:22].'+3)
+
+%!function r = f(c, x, y)
+%!  r = c{2} + x + y;
+%!endfunction
+%!
+%!test % setting IdxDimensions, with non-indexed cell array
+%! c = {1, 2};
+%! x = 1:5;
+%! y = 2:6;
+%! r = pararrayfun(2, @f, c, x, y, "VerboseLevel", 1, "IdxDimensions", [0 2 2]);
+%! assert(r, (5:2:13).')
+
+%!test % no output
+%! f = @(x) disp(x(:)');
+%! pararrayfun(2, f, [1:40], "Vectorized", true, "ChunksPerProc", 4, "VerboseLevel", 1)
+
+%!demo % see how it works
+%! f = @(x) disp(x(:)');
+%! pararrayfun(2, f, [1:40], "Vectorized", true, "ChunksPerProc", 4, "VerboseLevel", 2)
diff --git a/parcellfun.m b/parcellfun.m
--- a/parcellfun.m
+++ b/parcellfun.m
@@ -75,24 +75,50 @@
   endif
 
   [nargs, uniform_output, error_handler, ...
-  verbose_level, chunks_per_proc] = parcellfun_opts (varargin);
+    verbose_level, chunks_per_proc, ...
+    ~, cat_dimensions, idx_dimensions] = parcellfun_opts (varargin);
 
   args = varargin(1:nargs);
-  if (! all (cellfun ("isclass", args, "cell")))
-    error ("parcellfun: all non-option arguments except the first one must be cell arrays");
-  endif
+    
+  # use default indexing ?
+  default_indexing = isempty(idx_dimensions);
 
   if (nargs == 0)
     print_usage ();
-  elseif (nargs > 1)
-    [err, args{:}] = common_size (args{:});
-    if (err)
-      error ("parcellfun: arguments size must match");
+  elseif default_indexing
+    if (nargs > 1)
+      [err, args{:}] = common_size (args{:});
+      if (err)
+        error ("parcellfun: arguments size must match or\n\
+                   need IdxDimensions option");
+      endif
     endif
+    indexed_args_mask = ones(size(args), "logical");
+    njobs = numel (args{1});
+  else
+    # FIXME: code duplicated between parcellfun, pararrayfun, and chunk_parcellfun
+    assert(numel(idx_dimensions), numel(args));
+    # indexes of the indexed arguments
+    indexed_args_mask = (idx_dimensions != 0);
+    # sizes of the indexed arguments along the indexed dimension
+    indexed_dimensions_sizes = cellfun(@(a, dim) size(a, dim), 
+                                       args(indexed_args_mask), 
+                                       num2cell(abs(idx_dimensions(indexed_args_mask))),
+                                       "UniformOutput", false);
+    # all these lengths must be equal
+    if numel(indexed_dimensions_sizes) > 1
+      assert(isequal(indexed_dimensions_sizes{:}))
+    endif
+    # number of "jobs" = length along the indexed dimension
+    # they are all equal, just take the first one
+    njobs = indexed_dimensions_sizes{1};
   endif
 
-  njobs = numel (args{1});
-
+  if (! all (cellfun ("isclass", args(indexed_args_mask), "cell")))
+    error ("parcellfun: all non-option arguments except the first one must be cell arrays or\n\
+            need IdxDimensions option");
+  endif
+  
   if (chunks_per_proc > 0 && chunks_per_proc < njobs / nproc)
     ## We need chunked evaluation.
 
@@ -105,11 +131,12 @@
     endif
 
     [varargout{1:nargout}] = chunk_parcellfun (nproc, chunks_per_proc, ...
-    chunk_fun, [], verbose_level, args{:});
-    return
+      chunk_fun, [], verbose_level, cat_dimensions, idx_dimensions, args{:});
+    return ## return
   endif
 
-  nproc = min (nproc, numel (args{1}));
+  ## no need for all the procs if fewer jobs to do
+  nproc = min (nproc, njobs);
 
   ## create communication pipes.
   cmdr = cmdw = resr = resw = zeros (nproc, 1);
@@ -133,6 +160,9 @@
   if (err)
     error ("failed to open pipe: %s", msg);
   endif
+  
+  # backup njobs (will be changed later)
+  njobs_total = njobs;
 
   iproc = 0; # the parent process
   nsuc = 0; # number of processes succesfully forked.
@@ -226,9 +256,30 @@
             res = cell (1, nargout);
             argsc = cell (1, nargs);
             for i = 1:nargs
-              argsc{i} = args{i}{cmd};
+              if ( default_indexing )
+                argsc{i} = args{i}{cmd};
+              elseif (idx_dimensions(i) != 0)
+                # prepare indexing
+                idx = {":"}(ones(1, ndims(args{i})));
+                # keep only the values corresponding to the cmd index in the indexed dimension
+                idx(abs(idx_dimensions(i))) = cmd;
+                if (idx_dimensions(i) < 0)
+                  ## no reshape
+                  argsc{i} = [args{i}{idx{:}}];
+                else
+                  ## reshape
+                  # extract the values
+                  # get the template for the shape
+                  shape = size(args{i});
+                  # replace the indexed dim by 1
+                  shape(idx_dimensions(i)) = 1;
+                  argsc{i} = reshape([args{i}{idx{:}}], shape);
+                endif
+              else
+                # argument not indexed, pass it "as is"
+                argsc{i} = args{i};
+              endif
             endfor
-
             if (isempty (error_handler))
               ## unguarded evaluation.
               [res{:}] = fun (argsc{:});
@@ -259,7 +310,7 @@
 
         ## just indicate the error. don't quit this function !!!!
         fputs (stderr, "\n");
-        warning ("parcellfun: unhandled error in subprocess %d", iproc);
+        warning ("parcellfun: unhandled error in subprocess %d\n  %s", iproc, lasterr ());
 
         ## send a termination notice.
         fwrite (statw, -iproc, "double");
@@ -372,6 +423,7 @@
 
     end_unwind_protect
 
+    if false # FIXME to be removed
     ## we're finished. transform the result.
     varargout = cell (1, nargout);
     shape = size (varargin{1});
@@ -381,7 +433,100 @@
         varargout{i} = cell2mat (varargout{i});
       endif
     endfor
-
+    endif
+  
+    ## Concatenate output args and reshape them to correct size.
+    ## set the requested sizes of the output arguments
+    if isempty(cat_dimensions)
+      if default_indexing
+        shape = size (args{1});
+      else
+        ## output as column vector
+        shape = [njobs_total 1];
+      endif
+      ## duplicate shape
+      shapes = {shape}(ones(nargout, 1));
+      cat_dimensions = ones(1, nargout);
+      # # bring all outputs to column vectors
+      # res = cellfun(
+        # @(argout) cellfun(@(chunk_output) chunk_output(:), argout, ...
+                          # "UniformOutput", false),
+        # res, "UniformOutput", false);
+    else
+      # FIXME : (from pararrayfun) never tested, probably not working yet
+      assert(numel(cat_dimensions), nargout);
+      ## take each argout first chunk, 
+      ## return its size with the corresponding catdimension'th size replaced by 
+      ## the number of elements along the indexed dims
+      shapes = cell(nargout, 1);
+      for idx_argout = 1:nargout
+        # the results are splitted, take the first chunk as a template
+        shape = ones(1, cat_dimensions(idx_argout));
+        # replace the right dim of shape with njobs_total
+        shape(cat_dimensions(idx_argout)) = njobs_total;
+        shapes{idx_argout} = shape;
+        # in case of uniform_output, cell2mat will do the actual concatenation
+      endfor
+    endif  
+    ## for each argout concatenate the chunks along its cat_dimension, 
+    ## and then reshape to its shape
+    ## remember : res = cell (nargout, njobs_total);
+    ## can not loop directly over argouts like in chunk_parcellfun, but over indexes
+    varargout = cellfun (@(argout_idx, shape) reshape (res(argout_idx, :), shape), 
+      num2cell((1:nargout).'), shapes, "UniformOutput", false );
+      
+    if (uniform_output)
+      varargout = cellfun(@(c) cell2mat(c), varargout, "UniformOutput", false );
+    endif
   endif
 
 endfunction
+
+# tests borrowed from libinterp/corefcn/cellfun.cc
+
+%!assert(parcellfun(2, @(x) x, {1, 2}), [1, 2])
+%!assert(parcellfun(2, @(x) x, {1, 2}, 'uniformoutput', false), {1, 2})
+
+%% Number of input and output arguments may not be limited to one
+%!test
+%! A = parcellfun (2, @(x,y,z) x + y + z, {1, 1, 1}, {2, 2, 2}, {3, 4, 5});
+%! assert (A, [6, 7, 8]);
+%!test
+%! A = parcellfun (2, @(x,y,z) x + y + z, {1, 1, 1}, {2, 2, 2}, {3, 4, 5}, \
+%!              "UniformOutput", false);
+%! assert (A, {6, 7, 8});
+
+%!test %% Two input arguments of different types
+%! A = parcellfun (2, @(x,y) islogical (x) && ischar (y), {false, true}, {"a", 3});
+%! assert (A, [true, false]);
+
+%!test %% Pass another variable to the anonymous function
+%! y = true;
+%! A = parcellfun (2, @(x) islogical (x) && y, {false, 0.3});
+%! assert (A, [true, false]);
+
+%!test %% Three ouptut arguments of different type
+%! [A, B, C] = parcellfun (2, @find, {10, 11; 0, 12}, "UniformOutput", false);
+%! assert (isequal (A, {true, true; [], true}));
+%! assert (isequal (B, {true, true; [], true}));
+%! assert (isequal (C, {10, 11; [], 12}));
+
+
+# other tests
+
+%!test % UniformOutput = true
+%! f = @(x, y) x + y;
+%! x = [1:10; 2:11];
+%! m = num2cell(x);
+%! r = parcellfun(2, f, m, m, "UniformOutput", true);
+%! assert(r, f(x, x))
+%! r = parcellfun(2, f, m, m);
+%! assert(r, f(x, x))
+
+%!test % IdxDimensions
+%! f = @(x, y) sum(x + y.');
+%! x = [1:10; 2:11];
+%! m = num2cell(x);
+%! r = parcellfun(2, f, m, m.', "IdxDimensions", [2 1]);
+%! assert(r, f(x, x.') (:) )
+
diff --git a/private/chunk_parcellfun.m b/private/chunk_parcellfun.m
--- a/private/chunk_parcellfun.m
+++ b/private/chunk_parcellfun.m
@@ -20,34 +20,112 @@
 ## @end deftypefn
 
 function varargout = chunk_parcellfun (nproc, chunks_per_proc, func,
-  error_handler, verbose_level, varargin)
+  error_handler, verbose_level, cat_dimensions, idx_dimensions, varargin)
 
   args = varargin;
   
   nchunks = chunks_per_proc * nproc;
+  
+  # use default indexing ?
+  default_indexing = isempty(idx_dimensions);
 
   ## Compute optimal chunk sizes.
-  N = numel (args{1});
+  if default_indexing
+    idx_dimensions = ones(1, numel(args));
+    N = numel (args{1});
+    ## make all arguments column vectors, to prepare for mat2cell
+    args = cellfun(@(argument) argument(:), args, "UniformOutput", false);
+    # all indexed
+    indexed_args_mask = ones(size(args), "logical");
+  else
+    ## find the sizes of the input arguments 
+    ## along the corresponding indexing dimensions
+    indexed_args_mask = (idx_dimensions != 0);
+    indexed_dimensions_sizes = cellfun(@(array, dim) size(array, dim),  
+                             args(indexed_args_mask),
+                             num2cell(idx_dimensions(indexed_args_mask)));
+    ## check that they are all the same
+    assert(all(indexed_dimensions_sizes == indexed_dimensions_sizes(1)), 
+           "parcellfun: input arguments indexing dimension size mispatch");
+    N = indexed_dimensions_sizes(1);
+  endif
   len_chunk = ceil (N/nchunks);
   chunk_sizes = len_chunk (ones(1, nchunks));
   chunk_sizes(1:nchunks*len_chunk - N) -= 1;
 
+  ## set the split sizes
+  split_sizes = cellfun(@(argument) num2cell(size(argument)), 
+                          args, "UniformOutput", false);
+  for arg_idx = 1:numel(args)
+    if (idx_dimensions(arg_idx) != 0)
+      split_sizes{arg_idx}{idx_dimensions(arg_idx)} = chunk_sizes;
+    endif
+  endfor
   ## Split argument arrays into chunks (thus making arrays of arrays).
-  chunked_args = cellfun (@(arg) mat2cell (arg(:), chunk_sizes), args, ...
-  "UniformOutput", false);
+  chunked_args = cell(size(args));
+  # split only indexed arguments
+  # idx = find(indexed_args_mask)
+  chunked_args(indexed_args_mask) = cellfun (
+    @(argument, split_size) mat2cell (argument, split_size{:}), ...
+      args(indexed_args_mask), split_sizes(indexed_args_mask), ...
+      "UniformOutput", false);
+  # not indexed arguments are passed "as is"
+  chunked_args(! indexed_args_mask) = args(! indexed_args_mask);
+  # chunked_args = cellfun (
+    # @(argument, split_size) mat2cell (argument, split_size{:}), ...
+      # args, split_sizes, ...
+      # "UniformOutput", false);
 
   ## Attach error handler if present.
   if (! isempty (error_handler))
     chunked_args = [chunked_args, {"ErrorHandler", error_handler}];
   endif
-
+  # pass additional options
+  chunked_args = [chunked_args, {"UniformOutput", false, "VerboseLevel", verbose_level}];
+  # set IdxDimensions only if required
+  if (! default_indexing)
+    chunked_args = [chunked_args, {"IdxDimensions", -idx_dimensions}];
+  endif
   ## Main call.
-  [out_brut{1:nargout}] = parcellfun (nproc, func, chunked_args{:}, ...
-  "UniformOutput", false, "VerboseLevel", verbose_level);
+  ## use the original idx_dimensions (=> ensure the original behaviour)
+  [out_brut{1:nargout}] = parcellfun (nproc, func, chunked_args{:});
 
   ## Concatenate output args and reshape them to correct size.
-  true_size = size (args{1});
-  varargout = cellfun (@(arg) reshape (vertcat (arg{:}), true_size), out_brut, "UniformOutput", false);
-
+  ## set the requested sizes of the output arguments
+  if isempty(cat_dimensions)
+    if default_indexing
+      true_size = size (varargin{1});
+    else
+      ## output as column vector
+      true_size = [N 1];
+    endif
+    ## duplicate true_size
+    true_sizes = {true_size}(ones(1, nargout));
+    cat_dimensions = ones(1, nargout);
+    # bring all outputs to column vectors
+    out_brut = cellfun(
+      @(argout) cellfun(@(chunk_output) chunk_output(:), argout, ...
+                        "UniformOutput", false),
+      out_brut, "UniformOutput", false);
+  else
+    assert(numel(cat_dimensions), nargout);
+    ## take each argout first chunk, 
+    ## return its size with the corresponding catdimension'th size replaced by 
+    ## the number of elements along the indexed dims
+    true_sizes = cellfun(
+       @(argout, dim) [ size(argout{1, 1})(1:(dim-1)) , ...
+                        N , ...
+                        size(argout{1, 1})((dim+1):ndims(argout{1, 1})) ], 
+       out_brut, num2cell(cat_dimensions), "UniformOutput", false );
+  endif 
+  if ! isempty(out_brut)
+    ## for each argout concatenate the chunks along its cat_dimension, 
+    ## and then reshape to its true_size
+    varargout = cellfun (
+      @(arg, true_size, cat_dimension) 
+        reshape (cat(cat_dimension, arg{:}), true_size), 
+      out_brut, true_sizes, num2cell(cat_dimensions), "UniformOutput", false );
+  endif
+  
 endfunction
 
diff --git a/private/parcellfun_opts.m b/private/parcellfun_opts.m
--- a/private/parcellfun_opts.m
+++ b/private/parcellfun_opts.m
@@ -19,14 +19,17 @@
 ## @end deftypefn
 
 function [nargs, uniform_output, error_handler, ...
-  verbose_level, chunks_per_proc, vectorized] = parcellfun_opts (args)
+  verbose_level, chunks_per_proc, vectorized, ...
+  cat_dimensions, idx_dimensions] = parcellfun_opts (args)
 
   uniform_output = true;
   error_handler = [];
   verbose_level = 1; # default to normal output level
   chunks_per_proc = 0; # 0 means than size of chunk is 1
   vectorized = false;
-
+  cat_dimensions = [];
+  idx_dimensions = [];
+  
   nargs = length (args);
 
   ## parse options
@@ -64,8 +67,14 @@
       if (! isscalar (vectorized))
         error ("parcellfun: Vectorized must be a logical scalar");
       endif
+    case "catdimensions"
+      assert(isvector(val), "parcellfun: CatDimensions must be a vector");
+      cat_dimensions = val;
+    case "idxdimensions"
+      assert(isvector(val), "parcellfun: IdxDimensions must be a vector");
+      idx_dimensions = val;
     otherwise
-      break;
+      error(sprintf('parcellfun: unknown option "%s"', opt));
     endswitch
     nargs -= 2;
   endwhile
