Clear["`*"]

V  =  
      1/Sqrt[r^2 +  a^2  - 2 a r μx ]  -  
        1/Sqrt[r^2 +  a^2  +  2 a r μx ]   +   1/
          Sqrt[r^2 +  a^2  - 2 a r μy ]  -  1/Sqrt[
            r^2 +  a^2  +  2 a r μy ]

Series[V, {r, Infinity, 4}]

l = 1

Table[
      SphericalHarmonicY[l, m, π/2, 0] + 
        SphericalHarmonicY[l, m, π/2, π/2], {m, -l, 
        l}]

l  =  3

Do[t[l] =  
        Table[2  
            q a^l (SphericalHarmonicY[l, m, π/2, 0] + 
                SphericalHarmonicY[l, m, π/2, -π/2]), {m, 0, 
            l}], {l, 5}]

t[1]

t[3]

t[5]

Do[tt[l] =  
        Table[SphericalHarmonicY[l, m, 0, π/2], {m, 0, l}], {l, 1, 
        10}]

tt[1]

Table[tt[l], {l, 2, 10, 2}] // TableForm

Plot in x-y plane

Clear["`*"]

Φ  =  1/r^3 LegendreP[2, μ]

V[x_, 
        z_]  =  Φ  /.  {r ->  
            Sqrt[x^2 + z^2],  μ ->  z/Sqrt[x^2 + z^2]}

V[2, -1]

ContourPlot[V[x, z], {x, -2, 2}, {z, -2, 2}, 
      PlotPoints ->  200]

cp  = 
      ContourPlot[V[x, z], {x, -2, 2}, {z, -2, 2}, 
        PlotPoints -> 200, 
        DisplayFunction ->  Identity]

cir = Graphics[{Thickness[0.008], Circle[{0, 0}, 1]}]

Show[{cp, cir}, DisplayFunction ->  $DisplayFunction]

Clear["`*"]

V =  
      1/Sqrt[x^2 + (z - a)^2] +  1/Sqrt[x^2 + (z + a)^2]  - 
        2/Sqrt[x^2 + z^2]

V0  =  Normal[Series[V, {a, 0, 3}]]

Vp  =  V /.  {a -> 1, z ->  0}

V0p  =  V0 /.  {a ->  1,  z ->  0}

qq =  
      Plot[{V0p, Vp}, {x, 1, 3}, 
        PlotLabel ->  "Potential (Black) vs. Expansion (Red)", 
        AxesLabel ->  {"r", ""}, 
        PlotStyle ->  {{Thickness[0.01], RGBColor[1, 0, 0]}, 
            Thickness[0.01]}, 
PlotRange ->  {{0.5, 
              3}, {-1, 0}}]

Export["C:\Documents and Settings\johnson\My \
Documents\Filexfer\Class06S\prob41.eps", qq]