Finding possible arrangements of books on a shelf with Monte-Carlo in Matlab

I am trying to learn about Monte-Carlo simulations and data structures in Matlab, and as an example to play with I am trying to estimate numerically in Matlab how many possible ways books can be arranged on a shelf.

I have x3 categories "physics", "sci-fi", and "travel". Each contains N_phys, N_scifi, and N_travel numbers of books respectively. Within their category, the books can be placed in any order, but they must stay within their respective category (i.e all the physics books together). The categories can then be arranged in any order (i.e I could have all the sci-fi books first, followed by travel, followed by physics, for example).

This problem can be solved easily by a permutation, and the result is N_phys ! x N_scifi ! x N_travel ! x 3 !. For N_phys = 4, N_scifi = 3, and N_travel = 2, the result is 1728 possible arrangements.

I have decided to label both each book and its category by integers, with "physics" = 1, "sci-fi" = 2, and "travel" = 3. The shelf is then a 2D array. So for example, a whole shelf could look like the following:

3   2   1   4   1   2   2   1   3
1   1   1   1   3   3   2   2   2

where the the first 4 books are physics (because the second row is 1), followed by 2 travel books, and finally 3 sci-fi books, like this:

I have written a "brute-force" approach, shown below, which seems to give the correct result:

N_phys = 4;   % Number of physics books
N_scifi = 3;  % Number of sci-fi books
N_travel = 2; % Number of travel books

books_physics = [1:N_phys;...
                 ones(1,N_phys)]; % Collection of physics books

books_scifi = [1:N_scifi;...
               2*ones(1,N_scifi)]; % Collection of sci-fi books

books_travel = [1:N_travel;...
                3*ones(1,N_travel)];  % Collection of travel books

num_samples = 10000; % Number of random shelves to generate
unique_shelves = zeros(num_samples*2,N_phys+N_scifi+N_travel); % Preallocate 

unique_shelves(1:2,:) = [books_physics books_scifi books_travel];
num_unique_shelves = 1;

% Generate "num_samples" permutations of shelves randomly
for sample_num = 1:num_samples
    
books_physics_shuffled = books_physics( :, randperm(N_phys) ); % Shuffle physics books
books_scifi_shuffled = books_scifi( :, randperm(N_scifi) );    % Shuffle sci-fi books
books_travel_shuffled = books_travel( :, randperm(N_travel) ); % Shuffle travel books

category_order = randperm(3,3); % Choose order of categories, e.g. sci-fi/phsycis/travel  = [2 1 3]

shelf = [];

% Arrange the categories in a random order
for k = 1:3
    if category_order(k) == 1
        shelf = [shelf books_physics_shuffled];
    elseif category_order(k) == 2
        shelf = [shelf books_scifi_shuffled];
    elseif category_order(k) == 3
        shelf = [shelf books_travel_shuffled];
    end
end

% Iterate over discovered shelves, and see if we have found a new unique one
shelf_exists = 0;
for k = 1:num_unique_shelves
    if shelf == unique_shelves( (2*k-1):(2*k),:)
        shelf_exists = 1; % Shelf was previously discovered
        break
    end
end

if ~shelf_exists % New shelf was found
    unique_shelves( (2*num_unique_shelves+1):(2*num_unique_shelves+2),:) = shelf; % Add shelf to existing ones
    num_unique_shelves = num_unique_shelves + 1;
end

end

disp(['Number of unique shelves found = ',num2str(num_unique_shelves)])
disp(['Expected = ', num2str(factorial(N_phys)*factorial(N_scifi)*factorial(N_travel)*factorial(3) )])

As can be seen, I am basically randomly generating a shelf, and then checking if it has been previously found. If it has, I add it to the list.

My question is how can I learn and improve this approach, making it more concise? Is there a better data structure for storing such "unique_shelves", instead of the 2D array of integers, as well as finding the matches?

This is also quite slow to run, and also doesn't scale for more categories, since they are hardcoded.

Tips or alternative examples would be great! Thanks.