The acoustic black hole (ABH) effect can be achieved embedding cuneate indentations with power-law profile in plates. That results in remarkable properties like the reduction of plate vibrations, the potential for energy harvesting thanks to the focalization of energy, or the design of new metamaterials to manipulate acoustics waves. The analysis of such phenomena demands simulating the modal shapes and response to external excitations of ABH plates. This is usually done by means of numerical approaches, like the finite element method (FEM). However, if one is interested in performing long parametric analyses and in capturing the ABH plate behavior at high frequencies, the computational cost associated to numerical methods may become too demanding. In this work, a semi-analytical approach is suggested to circumvent the situation. The Rayleigh-Ritz method is applied using two-dimensional Gaussian functions to expand the flexural motion of a plate with non-uniform thickness. Then, a matrix-replacing strategy is proposed to embed the multiple ABHs in the plate. That results in low dimensional matrix systems, which yet provide very accurate solutions. After presenting all theoretical developments, the semi-analytical method is first applied to analyze the performance of a single ABH when varying several of its parameters. Then, various configurations involving multiple ABHs are considered. Those range from long strips exhibiting frequency attenuation bands, to plates containing ABH indentations in many shapes and sizes.