Solution (source code)

= Solution

For large $|z|$, expand in $a^2/z^2$ to obtain
$$
\boxed{\psi(0,z)=\frac{Q}{4\pi\epsilon_0}\left[\frac1{|z|}-\frac{a^2}{4|z|^3}+O(a^4/|z|^5)\right].}
$$
The axisymmetric exterior <multipole expansion> replaces each $|z|^{-n-1}$ on the positive axis by $r^{-n-1}P_n(\cos\theta)$. Reflection symmetry eliminates the odd multipoles, so
$$
\boxed{\psi(r,\theta)=\frac{Q}{4\pi\epsilon_0}\left[\frac1r-\frac{a^2}{4r^3}P_2(\cos\theta)+O(a^4/r^5)\right],\qquad r\gg a.}
$$
The leading term is the potential of total charge $Q$ concentrated at the center; the next is its <electric quadrupole> correction. On either axis $P_2(\pm1)=1$, recovering the same even axis potential.