A universal algorithm for filling an array

I have a code and it works for 3 dimensions, but I want to be able to specify the number of dimensions dynamically. There is a function nchoosek, but it creates an array of non-repeating elements.

This is my code for 3D:

energy_band=70;
count=1;
for i=-energy_band:energy_band
    for j=-energy_band:energy_band
        for z=-energy_band:energy_band
          out(1,count)=i;
          out(2,count)=j;
          out(3,count)=z;
          count=count+1;
       end
   end
end

We need an analogue of this algorithm, but for an arbitrary number of dimensions. For example, when I need to create an array for 4 dimensions, the code will look like this:

for i=-energy_band:energy_band %array formation cycle with energy shifts in the range from -1.75 eV to 1.75 eV
            for j=-energy_band:energy_band
                for z=-energy_band:energy_band
                    for k=-energy_band:energy_band
                        out(1,count)=i;
                        out(2,count)=j;
                        out(3,count)=z;
                        out(4,count)=k;
                        count=count+1;
                    end
                end
            end
        end

At the moment I am solving this problem with the following function, but as you can see there is no universality here:

function out=energy_array(matrix_version,energy_band)
count=1;
switch matrix_version
    case 1
        for i=-energy_band:energy_band %array formation cycle with energy shifts in the range from -1.75 eV to 1.75 eV
            out(1,count)=i;
            count=count+1;
        end
    case 2
        for i=-energy_band:energy_band %array formation cycle with energy shifts in the range from -1.75 eV to 1.75 eV
            for j=-energy_band:energy_band
                out(1,count)=i;
                out(2,count)=j;
                count=count+1;
            end
        end
    case 3
        for i=-energy_band:energy_band %array formation cycle with energy shifts in the range from -1.75 eV to 1.75 eV
            for j=-energy_band:energy_band
                for z=-energy_band:energy_band
                    out(1,count)=i;
                    out(2,count)=j;
                    out(3,count)=z;
                    count=count+1;
                end
            end
        end
    case 4
        for i=-energy_band:energy_band %array formation cycle with energy shifts in the range from -1.75 eV to 1.75 eV
            for j=-energy_band:energy_band
                for z=-energy_band:energy_band
                    for k=-energy_band:energy_band
                        out(1,count)=i;
                        out(2,count)=j;
                        out(3,count)=z;
                        out(4,count)=k;
                        count=count+1;
                    end
                end
            end
        end
    case 5
        for i=-energy_band:energy_band %array formation cycle with energy shifts in the range from -1.75 eV to 1.75 eV
            for j=-energy_band:energy_band
                for z=-energy_band:energy_band
                    for k=-energy_band:energy_band
                        for l=-energy_band:energy_band
                            out(1,count)=i;
                            out(2,count)=j;
                            out(3,count)=z;
                            out(4,count)=k;
                            out(5,count)=l;
                            count=count+1;
                        end
                    end
                end
            end
        end
end

How can I implement this more elegantly, without spelling out all the different cases?

You can create such arrays using ndgrid. This function outputs a separate array for each column of your desired array, so we need to apply some tricks to collect the arbitrary number of outputs and concatenate them as the columns of a matrix:

function out = energy_array(matrix_version, energy_band)
   arrays = cell(matrix_version,1);
   range = -energy_band:energy_band;
   inputs = repmat({range}, matrix_version, 1);
   [arrays{:}] = ndgrid(inputs{:});
   arrays = cellfun(@(x) x(:).', arrays, 'UniformOutput', false);
   out = cat(2, arrays{:});
end

arrays is a cell array with the number of outputs ndgrid will produce. This is the trick to receive a variable number of output arguments.

inputs is a cell array with the input arguments to ndgrid. inputs{:} is a comma-separated list of the values in the cell array, meaning each element is a separate argument.

arrays contains, after the call to ndgrid, a set of matrix_version arrays, each with matrix_version dimensions. We need to make these into column vectors, which you can do with arrays{i} = arrays{i}(:) for each i. This is what the cellfun line does. Finally, we concatenate the resulting arrays into a single matrix.

[Note: this is a better solution, I'm keeping this one around to prevent deletion.]