  We determine the full mass and $q^2$ dependence of the heavy quark vacuum
polarization function $\Pi(q^2)$ and its contribution to the total $e^+e^-$
cross section at ${\cal O}(\alpha_s^2)$ and ${\cal O}(\alpha_s^3)$ in
perturbative QCD. We use known results for the expansions of $\Pi(q^2)$ at high
energies, in the threshold region and around $q^2=0$, conformal mapping and the
Pad\'e approximation method. From our results for $\Pi(q^2)$ we determine
numerically at ${\cal O}(\alpha_s^3)$ the previously unknown non-logarithmic
contributions in the high-energy expansion at order $(m^2/q^2)^i$ for $i=0,1$
and the coefficients in the expansion around $q^2=0$ at order $q^{2n}$ with
$n\ge 2$. We also determine at ${\cal O}(\alpha_s^2)$ the previously unknown
${\cal O}(v^0)$ constant term in the expansion of $\Pi(q^2)$ in the threshold
region, where $v$ is the quark velocity. Our method allows for a quantitative
estimate of uncertainties and can be systematically improved once more
information in the three kinematic regions becomes available by future
multi-loop computations. For the contributions to the total $e^+e^-$ cross
section at ${\cal O}(\alpha_s^2)$ we confirm results obtained earlier by
Chetyrkin, K\"uhn and Steinhauser.
