We use the numerical renormalization group theory to investigate the Yu-Shiba-Rusinov (YSR) bound state properties of single magnetic molecules placed in an s-wave superconducting substrate. The molecule consists of a large core spin and a single orbital, coupled via exchange interaction. The critical Coulomb interaction for the singlet/doublet transition decreases in the presence of this exchange interaction for both Ferro and anti-ferromagnetic couplings. The number of YSR states also increases to two pairs, however, in the singlet phase, one of the pairs has zero spectral weight. To understand these results, we write down an effective model for the molecule in the limit of a large superconducting order parameter. Qualitatively, it explains the various phase transitions and spectral nature of the in-gap states.