We use variational methods to numerically study the matrix model for two-color QCD coupled to a single quark (matrix-QCD$_{2,1}$) in the ultra-strong Yang--Mills coupling ($g = \infty$). The model's spectrum is divided into superselection sectors characterized by baryon number ($B$) and spin ($J$). Our analysis focuses on the sectors ($B = 0, 1, 2$) and ($J = 0, 1$), which correspond to mesons, (anti-)diquarks, and (anti-)tetraquarks. For each sector, we examine the ground-state properties in both the chiral and heavy-quark limits, revealing a rich structure of quantum phase transitions (QPTs). We further analyze how the total spin is distributed between the gluonic and quark components, and find that the gluonic contribution is substantial in several sectors. Finally, by introducing a finite baryon chemical potential ($\mu$), we construct the model's phase diagram and demonstrate that, for sufficiently large ($\mu$), the ground state can acquire non-zero spin.