In this talk I will explain how to obtain the perturbative NNNLO heavy quarkonium spectrum for equal and different masses. This computation allows to determine the charm and bottom quark masses from the bottomonium, charmonium and $B_c$ systems. The use of the renormalon subtracted scheme, provides control over the divergence of the perturbative series due to the pole mass renormalon. On top of this, we also study an alternative computational scheme that treats the static potential exactly.
Finally, I will present a determination of $\alpha_s(M_z)$ based on a renormalon free combination of the heavy quarkonium systems.