Prob 1

Clear["`*"]

Owing to the azimuthal symmetry of the charge density, the general expansion of the potential can be witten 
      Φ(r,θ) = (e/(4  π (Subscript[ϵ, 0]))) (∑_(l = 0)^∞) (Subscript[P, l])[cos[θ]] (Subscript[v, l])(
          r)  ∫(Subscript[P, l])[
              cos[θ']] Abs[
                  (Subscript[Y, 1  m])]^2 dΩ' 
where the Hartree screening potential (Subscript[v, l ]) is  given by  (Subscript[v, l])\(TraditionalForm\`\((r)\)\  = ∫\_0\%∞\(\( r\^l\)\_ < \)\/\(\(r\^\(l + 1\)\)\_ > \)\ \(Abs[
                  P2p(r)]\^2\) \(\(dr\)\(\ \)\(.\)\)\)
As a first step, we evaluate the angular integral. It is nonzero only for l=0 and l=2.

Ang[l_, m_]  := 
      2 π Integrate[
          Sin[θ] LegendreP[l, 
              Cos[θ]] SphericalHarmonicY[1, m, θ, 
                0]^2, {θ, 0, π}]

Table[Ang[0, m], {m, -1, 1}]

Table[Ang[2, m], {m, -1, 1}]

Now define the Hartree functions:

P2p[
        x_]  =  (1/(2 Sqrt[6])) ζ^(5/
              2) x^2 Exp[-ζ x/2]

P2p[
        x_]  =  (1/(2 Sqrt[6])) ζ^(5/
              2) x^2 Exp[-ζ x/2]

v[l_, r_]  :=  
      Simplify[(1/r^(l + 1)) 
            Integrate[x^l P2p[x]^2, {x, 0, r}, 
              Assumptions ->  ζ > 0]  + 
          r^l Integrate[P2p[x]^2/x^(l + 1), {x, r, Infinity}, 
              Assumptions ->  ζ > 0]]

v[0, r]

v[2, r]

Now we have the potential in general

Φ[m_, r_]  :=  (1/(4 π (Subscript[ϵ, 0]))) 
        Sum[ Ang[l, m] v[l, r] (Subscript[P, l])[μ], {l, 0, 2, 2}]

Φ[0, r]

Φ[1, r]

Table[
      Normal[Series[Φ[m, r], {r, 0, 3}]], {m, -1, 
        1}]

Table[
      Limit[Φ[m, r], ζ ->  Infinity, 
        Assumptions ->  r > 0], {m, -1, 1}]

Coefficient[v, Exp[-r ζ]]

v0[r_]  =  v0[r][[1]] +  
        Exp[- r ζ] Coefficient[v, 
            Exp[-r ζ]]

v2[r_]  =  
      Simplify[(1/r^3) 
            Integrate[x^2  P2p[x]^2, {x, 0, r}, 
              Assumptions ->  ζ > 0]  + 
          r^2 Integrate[P2p[x]^2/x^3, {x, r, Infinity}, 
              Assumptions ->  ζ > 0]]

v2[r_]  =  Expand[v2[r]]

Series[v2[r], {r, 0, 4}]

v  =  v2[r]  -  v2[r][[2]]

v2[r_]  =  v2[r][[2]] +  
        Exp[- r ζ] Coefficient[v, 
            Exp[-r ζ]]

Prob 2

A charge q is on the z axis at a distance d < a  from the center of a dielectric shell (inner radius a, outer radius b, dielectric constant ϵ/(Subscript[ϵ, 0]) = (Subscript[ϵ, r])) . 
(a)  Show that the leading terms in the asymptotic expansion of the potential in the region  r>b  are (1/(4  π (Subscript[ϵ, 0])))((q/r)+\(\(9\ \ d\ \ ϵ\)\/\(\(-2\)\ \((a/b)\)\^3\ \((\(-1\) + ϵ)\)\^2 + \ \((2 + 5\ ϵ + 2\ ϵ\^2)\)\)\) cos[θ]/(r^2)) . 
(b)  Find the electric field (direction and magnitude) at the origin. 
(c) Verify the correctness of your in the limit (Subscript[ϵ, r]) = 1.
(d) Verify the correctness of your in the special case b = a.

[In the following,we omit the coefficientq/4  π(Subscript[ϵ, 0])from the expansion of the point charge potential; this factor must be reintroduced into all final formulas .   We also designate(Subscript[ϵ, r])by ϵ.]

Clear["`*"]

Part a

Set up and solve equations for coefficients of (Subscript[P, 0])(cosθ)  and (Subscript[P, 1])(cosθ).  We designate the four equations for the coefficients of (Subscript[P, k])(cosθ) at the two boundaries by eq1(k), eq2(k), eq2(k), and eq4(k).

eq1[k_]   =  
      d^k/a^(2  k + 
                  1)  +  α  ==  β  +  γ/a^(2 k  + 1)

eq2[
        k_]  =  (k + 1) d^k/a^(2  k + 1)  - 
          k α  ==  ϵ (- 
                  k β  +  (k + 1) γ/
                  a^(2  k + 1) )

eq3[
        k_]  =  β  +  γ/
            b^(2 k + 1)  ==  δ/
          b^(2 k  +  1)

eq4[
        k_]  =  ϵ ( -k β  +  (k + 1) γ/
                  b^(2  k + 1) )  ==  (k + 1) δ/
            b^(2  k + 1)

eqs[k_]  =  {eq1[k], eq2[k], eq3[k], eq4[k]}

vars  =  {α, β, γ, δ}

sol  =  Solve[eqs[0], vars] // First

a[0]  =  FullSimplify[α /.  sol]

b[0]  =  FullSimplify[β /.  sol]

c[0]  =  FullSimplify[γ /.  sol]

d[
        0] = FullSimplify[δ /.  
          sol]

sol  =  Solve[eqs[1], vars] // First

a[1]  =  FullSimplify[α /.  sol]

b[1]  =  FullSimplify[β /.  sol]

c[1]  =  FullSimplify[γ /.  sol]

d[1]  =  FullSimplify[δ /.  sol]

Here are the first two terms in the expansion for r>b.  (This leads to the answer for part a)

ϕout   =  
      d[0]/r  +  Simplify[d[1]] Cos[θ]/r^2

Part b

Here are the first few terms (plus source) for the potential in the region r<a.

ϕin  =  
      1/Sqrt[r^2 + d^2 - 2 r d Cos[θ]]  +  a[0]  +  
        a[1] r Cos[θ]

Er  =  -Simplify[ D[ϕin, r] /.  r ->  0]

Eθ  =  - Simplify[D[ϕin, θ]/r]  /.  
        r -> 0

Here are the rectangular coordinates of the field at the origin

Ez  =  
      FullSimplify[Er Cos[θ]  - Eθ Sin[θ], 
        Assumptions ->  d > 0]

Ez = Simplify[Ez[[1]]] + Simplify[Ez[[2]]]

Ex  =  
      FullSimplify[
        Er Sin[θ]  +  Eθ Cos[θ]]

Part c

When ϵ = 1, the potential in the region r>a  reduces to the first two terms in the expansion of the potential of the point charge, and the field at the origin is just the fieldl of the point charge.

ϕout  /.  ϵ ->  1

Simplify[Ez, Assumptions -> d > 0] /.  ϵ ->  
        1

Part d

Again,when b=a, the potential in the region r>a is just the first two terms in the expansion of the potential of the point charge, and the field at the origin is just the fieldl of the point charge.

Expand[
      FullSimplify[ϕout, Assumptions  ->  { b = a}]]

Simplify[Ez, Assumptions -> d > 0] /.  b ->  a

In both cases  ϕout gives the first two terms in the asymptotic expansion
of the point charge potential, and Ez gives the field of the point charge at the origin.

Prob 3

A long cylindrical tube (inner radius a, outer radius b, relative permeability μ/(Subscript[μ, 0]))  with axis along z is placed in a uniform field of magnitude (Subscript[B, 0]) directed along the x axis. 
a) Determine the B field everywhere.
b) Show that the field inside the tube is uniform and find the factor by which it is reduced
 assuming that  μ/(Subscript[μ, 0]) = 5000 (99.8\% pure iron) and  b/a=5/4.
 c) Sketch the B-field.

Clear["`*"]

sol  =  
      Solve[{α  ==  β  + γ/
                  a^2, -α ==  μ (-β + γ/a^2),  β +  γ/
                  b^2  ==  -(Subscript[H, 0]) +  δ/
                  b^2,  μ ( -β  +  γ/b^2)  ==  
              (Subscript[H, 0])  +  δ/
                  b^2}, {α, β, γ, δ}] // 
        First

a1  =  FullSimplify[α /.  sol]

b1  =  FullSimplify[β /.  sol]

c1  =  FullSimplify[γ /.  sol]

d1  =  FullSimplify[δ /.  sol]

ϕin  =  a1 ρ Cos[ϕ]

ϕmid  =  (b1 ρ  +  
            c1/ρ)  Cos[ϕ]

ϕout  = ( -(Subscript[H, 0 ]) ρ   +  
            d1/ρ) Cos[ϕ]

Hinr  =  -D[ϕin, ρ]

Hinϕ  =  -D[ϕin, ϕ]/ρ

Hinx  =  
      Simplify[ Hinr Cos[ϕ] - Hinϕ Sin[ϕ]]

Hiny  =  
      Simplify[ Hinr Sin[ϕ] + Hinϕ Cos[ϕ]]

Hmidr  =  -D[ϕmid, ρ]

Hmidϕ  =  -D[ϕmid, ϕ]/ρ

Hmidx  =   Simplify[
        Hmidr Cos[ϕ] - Hmidϕ Sin[ϕ]]

dH2  =  
      Simplify[Coefficient[Hmidx, ρ^(-2)]]/ρ^2

Hmidx  =  Simplify[Hmidx - dH2]  +  dH2

Hmidy  =  
      Simplify[ Hmidr Sin[ϕ] + Hmidϕ Cos[ϕ]]

Houtr  =  -D[ϕout, ρ]

Houtϕ  =  -D[ϕout, ϕ]/ρ

Houtx  =   FullSimplify[
        Houtr Cos[ϕ] - Houtϕ Sin[ϕ]]

Houty =   FullSimplify[
        Houtr Sin[ϕ] +  Houtϕ Cos[ϕ]]

Simplify[Hinx/(Subscript[H, 0])  /.  b  ->  2  a]  /.  μ ->  
          5000 // N

Bx  =  
      Piecewise[{{Hinx, 0 < ρ < a}, {μ Hmidx, 
            a < ρ < b}, {Houtx, b < ρ}}]

By =  
      Piecewise[{{Hiny, 0 < ρ < a}, {μ Hmidy, 
            a < ρ < b}, {Houty, b < ρ}}]

<< Graphics`PlotField`

bx[x_, y_]  =  
      FullSimplify[
        Bx  /.  {ρ ->  Sqrt[x^2 + y^2], ϕ ->  
              ArcCos[y/Sqrt[x^2 + y^2]],  
            (Subscript[H, 0]) ->  1, 
a ->  1, 
            b ->  5/4}]

Clear[x, y]

by[x_, y_]  =  
      FullSimplify[
        By  /.  {ρ ->  Sqrt[x^2 + y^2], ϕ ->  
              ArcTan[y/x],  
            (Subscript[H, 0]) ->  1, 
a ->  1, 
            b ->  5/4}]

ff  =  
      PlotVectorField[{bx[x, y], by[x, y]} /.  μ ->  
            1000, {x, -2, 2}, {y, -2, 2}, PlotPoints ->  30, 
        ColorFunction ->  Hue, 
        DisplayFunction ->  Identity]

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

c2  =  
      Graphics[{Thickness[0.008], 
          Circle[{0, 0}, 5/4]}]

Show[{ff, c1, c2}, 
      DisplayFunction ->  $DisplayFunction]