If the variables are all the same 2D size, then the most efficient thing to do is load them one at a time and store the value into a slice of an array. Then once you have all of the data loaded, you can mean() across the third dimension.
dinfo = dir('matrix*.mat');
nfile = length(dinfo);
if nfile == 0
avg = [];
warning('There were no input files');
else
filestruct = load(dinfo(1).name);
fn = fieldnames(filestruct);
thisvar = filestruct.(fn{1});
orig_size = size(thisvar);
input_data = [thisvar(:), zeros(numel(thisvar), nfile-1, class(thisvar))];
for K = 2 : nfile
filestruct = load(dinfo(K).name);
fn = fieldnames(filestruct);
thisvar = filestruct.(fn{1});
if ~isequal(size(thisvar), orig_size)
error('file "%s" has a different size of matrix than the first file "%s"', dinfo(K).name, dinfo(1).name);
end
input_data(:,K) = thisvar(:);
end
input_data = reshape(input_data, [orig_size, nfile]);
avg = mean(input_data, ndim(input_data));
end
(with data of an unknown number of dimensions, it turned out to be easiest to store the data as columns and reshape afterwards.)
