In this research, we model Hulthén plus generalized inverse quadratic Yukawa potential to interact in a quark-antiquark system. The solutions of the Schrödinger equation are obtained using the Nikiforov-Uvarov method. The energy spectrum and normalized wave function were obtained. The masses of the heavy mesons for different quantum states such as 1S, 2S , 1P, 2P 3S, 4S, 1D, and 2D were predicted as 3.096 GeV, 3.686 GeV, 3.327 GeV, 3.774GeV, 4.040 GeV, 4.364GeV, 3.761 GeV, and 4.058 GeV respectively for charmonium (cc). Also, for bottomonium (bb) we obtained 9.460 GeV, 10.023 GeV, 9.841 GeV, 10.160 GeV, 10.345 GeV, 10.522 GeV, and 10.142GeV for different states of 1S , 2S , 1P , 2P , 3S , 4S , 1D respectively. The partition function was calculated from the energy spectrum, thereafter other thermal properties were obtained. The results obtained showed an improvement when compared with the work of other researchers and excellently agreed with experimental data with a percentage error of 1.60 % and 0.46 % for (cc) and (bb), respectively.