We determine the charm and bottom quark masses using the N 3 LO perturbative expression of the ground state (pseudoscalar) energy of the bottomonium, charmonium and B c systems: the η b , η c and B c masses. We work in the renormalon subtracted scheme, which allows us to control the divergence of the perturbation series due to the pole mass renormalon. Our result for the MS masses reads m c (m c ) = 1223(33) MeV and m b (m b ) = 4186(37) MeV. We also extract a value of α s from a renormalon-free combination of the η b , η c and B c masses: α s (M z ) = 0.1195(53). We explore the applicability of the weak coupling approximation to bottomonium n = 2 states. Finally, we consider an alternative computational scheme that treats the static potential exactly and study its convergence properties.