>> mygrad = @(m)gradient(permute(m,[2,1,3:ndims(m)]));
>> [gxx,gyx,gzx] = mygrad(mshx);
>> [gxy,gyy,gzy] = mygrad(mshy);
>> [gxz,gyz,gzz] = mygrad(mshz);
Note that this works correctly for all N, even the edge-case when N==1, i.e. a vector input:
mshx = ndgrid(0:.1:10);
Or you could just copy gradient Mfile to your working directory, save it with a different function name (e.g. mygrad), and remove the code that permutes the first two dimensions. It wouldn't be very hard.
"Why does gradient function calculate wrong dimensions 1 and 2?"
Because many users think of array horizontal dimension as x, and array vertical dimension as y... but these are the second and first dimensions respectively (of standard matrix/array size notation). This is probably best illustrated in the meshgrid help.
The gradient documentation does state: "The first output FX is always the gradient along the 2nd dimension of F, going across columns. The second output FY is always the gradient along the 1st dimension of F, going across rows. For the third output FZ and the outputs that follow, the Nth output is the gradient along the Nth dimension of F."
