"Problem 3.1"

Clear[a,b,V,Q]
" for even L ";

eq1=V Q[L]==2/(2 L+1)(A (a^L)+B (a^(-L-1)))
eq2=V Q[L]==2/(2 L+1)(A (b^L)+B (b^(-L-1)))

FullSimplify[Flatten[Solve[{eq1,eq2},{A,B}]]]
Aeven[L_]=A /. %
Beven[L_]=B /. %%

" for odd L ";

eq1=V Q[L]==2/(2 L+1)(A (a^L)+B (a^(-L-1)))
eq2=-V Q[L]==2/(2 L+1)(A (b^L)+B (b^(-L-1)))

FullSimplify[Flatten[Solve[{eq1,eq2},{A,B}]]]
Aodd[L_]=A /. %
Bodd[L_]=B /. %%

Table[Integrate[LegendreP[n,x],{x,0,1}],{n,0,10}]

" from Eq. (3.26) in text; only for odd L";
Q[L_]=(((-1/2))^((L-1)/2)) (L-2)!!/2/(((L+1)/2)!)
" for L=0 ";
Q[0]=1;

" so for L=0 ";

Aeven[0]
Beven[0]

Table[FullSimplify[Aodd[k]],{k,1,3,2}]
Table[FullSimplify[Bodd[k]],{k,1,3,2}]

"project onto z-x plane";
r=Sqrt[x^2+y^2+z^2] /. y->0
Cosθ=z/Sqrt[x^2+y^2+z^2]  /. y->0

V=1;
a=1;
b=2;
Lmax=21;

phi[x_,z_]=N[Aeven[0] LegendreP[0,Cosθ]+Sum[(Aodd[L] (r^L) +Bodd[L] 1/(r^(L+1 ))) LegendreP[L,Cosθ],{L,1,Lmax,2}]];

ContourPlot[phi[x,z],{x,-2,2},{z,-2,2},RegionFunction->Function[{x,z},a<=(Sqrt[x^2+z^2])≤b  ],Contours->{0,0.1,0.25,0.5,1},ContourStyle->{{Green,Thickness[0.01]},Purple,Cyan,Yellow,{Blue,Thickness[0.01]}},ContourShading->None]

Plot3D[phi[x,z],{x,-2,2},{z,-2,2},RegionFunction->Function[{x,z},a<(Sqrt[x^2+z^2])<b]]