Comparison of potential obtained by image charges to the
the potential obtained by expansion in spherical harmonics.

F1[r_,th_, a_,lx_]:=Sum[r^l/a^(l+1)*Sqrt[4*Pi/(2*l+1)]*SphericalHarmonicY[l,0,th,0],{l,1,lx,2}];
F2[r_,th_, a_,lx_]:=Sum[a^l/r^(l+1)*Sqrt[4*Pi/(2*l+1)]*SphericalHarmonicY[l,0,th,0],{l,1,lx,2}];
vim[r_,th_,a_]:=1/Sqrt[r^2+a^2-2*r*a*Cos[th]]-1/Sqrt[r^2+a^2+2*r*a*Cos[th]];

Note that the sum is only over the odd spherical harmonics

Dynamic@N[r]
th=Pi/3; lx=100; a=3/2;
t1=Table[{r,2F1[r,th,a,lx]},{r,1/40, 3,1/20}];
t2=Table[{r,2F2[r,th,a,lx]},{r,1/40, 3,1/20}];
tim=Table[{r,vim[r,th,a]},{r,1/40, 3,1/20}];
ListLinePlot[{t1,t2,tim},PlotRange->{{0,3},{-1,1}}]