We use the numerical renormalization group theory to investigate the Yu-Shiba-Rusinov (YSR) bound state properties of single magnetic molecules placed on an s-wave superconducting substrate. The molecule is modeled as 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 ferroand antiferromagnetic couplings. The number of YSR states also increases to two pairs; however, in the singlet phase, one of the pairs has zero spectral weight. We explore the evolution of the in-gap states using the Anderson model. Away from the particle-hole symmetry point, the results suggest a doublet-singlet-doublet transition as the on-site energy is lowered while keeping the Coulomb interaction fixed. We construct an effective model for the molecule to understand these results, in the limit of the large superconducting order parameter. Qualitatively, the model accounts for the phase transitions and spectral nature of the in-gap states. Finally, we analyze the effects of magnetic anisotropic fields of the core spin on in-gap states. Due to internal degrees of freedom of the spin excited states, a multitude of new states emerges within the gap. Depending on the sign and strength of the uniaxial anisotropic field, the results indicate up to three pairs of YSR states.