Clear["`*"]

In region 1:    ρ > b  (omitting the common factor Cos[θ])

ϕ1  =  -E0 ρ  +  c1/ρ

In region 2:    b > ρ > a

ϕ2  =  c2 ρ  +  c3/ρ

In region 3:  a  > ρ

ϕ3  =  c4 ρ

Determine the potential difference at a and b

Δϕa  =  ϕ2  -  ϕ3  /.  ρ ->  a

Δϕb  =  ϕ1  -  ϕ2  /.  ρ ->  b

Determine the discontinuity in D at the two interbaces

ΔDa  =  -D[ϕ3, ρ]  +  ϵ D[ϕ2, ρ]   /.  ρ ->  a

ΔDb  =  -ϵ D[ϕ2, ρ]  +   D[ϕ1, ρ]   /.  ρ ->  b

Set up and solve the boundary condition equations

eqns = {Δϕa == 
          0, Δϕb == 0, ΔDa == 
          0, ΔDb == 0}

sol  =  Solve[eqns, {c1, c2, c3, c4}] // First

c1= Simplify[c1 /.  sol]

c2 =  Simplify[c2 /.  sol]

c3  =  Simplify[c3 /.  sol]

c4  =  Simplify[c4 /.  sol]

Check that the solutions make sense

{c1, c2, c3, c4}  /.  ϵ ->  1

c1 // TeXForm

c2 //  TeXForm

c3  //  TeXForm

c4  //  TeXForm

Prepare a sketch for case ϵ = 2,  a=1,  b= 2,  E0=1

ϕs  =  
      Piecewise[{{ϕ1 Cos[θ], ρ > 
                2}, {ϕ2 Cos[θ], 
              2 > ρ > 1}, {ϕ3 Cos[θ],  
              1 > ρ > 0}}] /.  {ϵ -> 10, a -> 1, 
          b -> 2, E0 -> 1}

<< Graphics`ParametricPlot3D`

CylindricalPlot3D[ϕs, {ρ, 0, 2.5}, {θ, 0, 2  Pi}, 
      PlotLabel ->  "Dielectric Cylinder in External Field", 
      TextStyle -> {FontSize -> 12}]

ϕp  =  ϕs  /.  {ρ ->  
            Sqrt[y^2 + x^2], θ ->  
            ArcCos[x/Sqrt[x^2 + y^2]]}

cp  =  
      ContourPlot[ϕp, {x, -3, 3}, {y, -3, 3}, 
        PlotPoints ->  200, 
        DisplayFunction ->  Identity]

Needs["Graphics`PlotField`"]

ce  =  
      PlotGradientField[-ϕp, {x, -3, 3}, {y, -3, 3}, 
        ScaleFunction ->  (1 &),  
        DisplayFunction ->  Identity]

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

cir2 = Graphics[{Thickness[0.008], Circle[{0, 0}, 2]}]

Show[{cp, ce, cir1, cir2}, 
      DisplayFunction ->  $DisplayFunction]

Export["C:\Documents and Settings\johnson\My \
Documents\Filexfer\Class06S\prob48.eps", %66]