In this work, we investigate the quasibound states of charged massive scalar fields in the Kerr-Newman black hole spacetime by using a new method recently developed, which uses the polynomial conditions of the Heun functions. We calculate the resonant frequencies related to the spectrum of quasibound states and the wave functions, and also we analyze the (in)stability of the system. These results are particularized to the cases of the Schwarzschild and Kerr black holes. Additionally, we compare our analytical results with the numerical ones known in the literature. Finally, we apply this method to the supermassive black hole situated at the center of the M87 galaxy, compute the characteristic times of growth and decay of bosonic particles around the astrophysical object and discuss the results.