  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 $\eta_b$, $\eta_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{\rm MS}$ masses reads ${\overline m}_{c}({\overline
m}_{c})=1223(33)$ MeV and ${\overline m}_{b}({\overline m}_{b})=4186(37)$ MeV.
We also extract a value of $\alpha_s$ from a renormalon-free combination of the
$\eta_b$, $\eta_c$ and $B_c$ masses: $\alpha_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.
