If you have your matrices in a cell array instead using indices hidden in the name of the variables, the solution would be easier. But it works:
M = cat(4, Matrix1 + Matrix2 + Matrix3 + Matrix4);
M(isnan(M)) = 0;
meanM = sum(M, 4) ./ sum(M ~= 0, 4);
Instead of dividing the sum by the number of matrices, you divide by the number of non-zero elements. The NaNs vanish in the sum, if you replace them by zeros.
You had the right idea already, but stopped too early:
It's not possible to concatenate the 3D matrices
(unlike a 2D matrix) and calculate its mean.
