Solution (source code)

= Solution

The separated axisymmetric solution is
$$
\Phi(r,x)=\sum_{\ell=0}^{\infty}(A_\ell r^\ell+B_\ell r^{-\ell-1})P_\ell(x).
$$
Since $x(1-x^2)=\tfrac25(P_1-P_3)$, regularity at the origin and the <boundary condition> give
$$
\boxed{\Phi(r,x)=\frac25\left[\frac rR P_1(x)-\left(\frac rR\right)^3P_3(x)\right].}
$$