Generalising number of factorial input elements of a function

Dear Matlab community,

I would like to plot data representing subject-wise scores as beeswarm plots (showing variability across subjects + central tendency e.g. median), while labelling nested levels (groups&subgroups) into which the data is organised.

From various bits of code, I've compiled the function below, which is suited to a 2x2 division of the input data, as the example & screenshot below (of scores from 2 pieces, on 2 repetitions) illustrate.

Y = [15,17,5,7;13,15,3,5;16,18,6,8;14,16,4,6;13,15,3,5;17,19,7,9;15,17,5,7;13,15,3,5;15,17,5,7;17,19,7,9;14,16,4,6;15,17,5,7;19,21,9,11;12,14,2,4;16,18,6,8;14,16,4,6;17,19,7,9;15,17,5,7];
barAndBee (Y, 2, 2, 'median')

I would like to generalise this to any number of groups & subgroups in the data (Y), not just 2 groups and 2 subgroups. I know this task must somehow involve varargin, however given my limited knowledge of Matlab, adapting the code is proving very challenging. I am wondering if someone can give me a few hints as to how I might go about it.

Full code of the barAndBee function is at https://pastebin.com/7TQAsTpd

Thanks!!