Electron-positron vacuum polarization potentials, accurate to all orders $n$ in ${(Z\ensuremath{\alpha})}^{n}$, are calculated using realistic nuclear charge distributions for selected nuclei with $26\ensuremath{\le}Z\ensuremath{\le}114$, and for the quasimolecule U-U. We give further details of the numerical methods used in calculations previously reported, and of the improvements which have led to the present increased accuracy. The vacuum polarization potentials are employed to determine muonic energy-level shifts [order ${(Z \ensuremath{\alpha})}^{n\ensuremath{\ge}1}\ensuremath{\alpha}$] for most relevant states with $l\ensuremath{\le}6$. For the shifts of order ${(Z \ensuremath{\alpha})}^{n\ensuremath{\ge}3}\ensuremath{\alpha}$, interpolation procedures for all values of $Z<137$ and for other exotic particles are noted. Although the effect of finite nuclear size is small for high-lying states (\ensuremath{\sim}10% reduction of the ${(Z \ensuremath{\alpha})}^{n\ensuremath{\ge}3}\ensuremath{\alpha}$ terms in heavy atoms), it yields a reduction by a factor \ensuremath{\gtrsim} 2 for low-lying states. In the case of the quasimolecule U-U, we find that the total vacuum polarization is a minor correction to the dominant Coulomb interaction, with the first-order (Uehling) contribution remaining much larger than the higher-order terms. We also present further details of a previous calculation of order ${(Z \ensuremath{\alpha})}^{2}{\ensuremath{\alpha}}^{2}$ in muonic lead.