We determine the charm and bottom quark masses using the N3LO 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 $\overline{MS}$ $$ \overline{\mathrm{MS}} $$ masses reads ${\overline{m}}_c\left({\overline{m}}_c\right)=1223\left(33\right)$ $$ {\overline{m}}_c\left({\overline{m}}_c\right)=1223(33) $$ MeV and ${\overline{m}}_b\left({\overline{m}}_b\right)=4186\left(37\right)$ $$ {\overline{m}}_b\left({\overline{m}}_b\right)=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.