Clear["Global`*"]
(* solve a cylindrical Poisson problem: 

a long rod is a dielectic with volume charge. in the infinite-length approximation the Laplace operator involves only s and phi, the z-derivatives of the potential vanish *)

(* © Marko Horbatsch, York University, Toronto, Canada *)

Q=2Pi Integrate[s Exp[-s^2/a^2],{s,0,Infinity}]

Poisson=D[V[x,y],{x,2}]+D[V[x,y],{y,2}]==-4Pi a^2/2 Exp[-(x^2+y^2)/a^2]

(* we set Q/L = λ =1 total linear charge density to 1 on the RHS *)

DSolve[Poisson,V,{x,y}]

(* an analytic solution in Cartesians doesn't work, so switch to cylindrical coordinates in which the equation is an ODE: *)

PoissonC=D[s D[Vc[s],s],s]/s==-4Pi a^2/2 Exp[-s^2/a^2]

solC=DSolve[PoissonC,Vc[s],s]

(* we need to fix two boundary conditions, and note that one is trivial since it is an additive constant C[1], but we cannot dispose of the other, since the log-divergence with s will give us the 1/s behaviour of the electric field at large distances. From physics we would argue that we want a finite potential at s=0: C[1]=-Pi a^4 *)

Limit[Vc[s]/.solC,s->0]

Vs=(Vc[s]/.solC/.{C[2]->0,C[1]->-Pi a^4}/.a->1/2) [[1]]

Plot[Vs,{s,0,10},PlotRange->All]

LogLinearPlot[Vs,{s,0.1,10},PlotRange->All]

Es=-D[Vs,s]

Plot[Es,{s,0,10},PlotRange->All]

(* a textbook solution that comes close: a piecewise constant density from s=0 to s=a, the field will rise linearly up to a peak value at s=a, and connect there to a 1/s fall-off *)

(* back to Cartesians: let us get a numerical solution using finite elements: we implement the Poisson equation on a rectangular (square) domain and apply what we know about the expected solution on the boundary *)

a=1/2

x0=5;
y0=5;

sol2d=NDSolveValue[{Poisson,V[x,y0]==-Pi a^4Log[Sqrt[x^2+y0^2]],V[x,-y0]==-Pi a^4Log[Sqrt[x^2+y0^2]],V[x0,y]==-Pi a^4Log[Sqrt[x0^2+y^2]],V[-x0,y]==-Pi a^4Log[Sqrt[x0^2+y^2]]},V,{x,-x0,x0},{y,-y0,y0}]

ContourPlot[sol2d[x,y],{x,-x0,x0},{y,-y0,y0},ContourLabels->All,Contours->{-0.3,-0.2,-0.1,0,0.1,0.15},Frame->True,FrameLabel->{x,y}]

(* Note: there are no deformations due to the use of Cartesian coordinates, the boundary condition was implemented on the boundary lines, and it is an appropriate BC. Now generalize the approach to use a disk geometry: *)

Ω=ImplicitRegion[((x)^2+(y)^2≤(10a)^2),{{x,-10a,10a},{y,-10a,10a}}];

(Subscript[Γ, D])=DirichletCondition[V[x,y]==-Pi a^4Log[10a],x^2+y^2==(10a)^2];

solI=NDSolveValue[{Poisson,(Subscript[Γ, D])},V,{x,y}∈Ω]

ContourPlot[solI[x,y],{x,-x0,x0},{y,-y0,y0},ContourLabels->All,Contours->{-0.3,-0.2,-0.1,0,0.1,0.15},Frame->True,FrameLabel->{x,y}]

(* how different is the solution when the BC is chosen as a zero reference value ? *)

(Subscript[Γ, D])=DirichletCondition[V[x,y]==0,x^2+y^2==(10a)^2];

solI0=NDSolveValue[{Poisson,(Subscript[Γ, D])},V,{x,y}∈Ω]

ContourPlot[solI0[x,y],{x,-x0,x0},{y,-y0,y0},ContourLabels->All(*,Contours->{-0.3,-0.2,-0.1,0,0.1,0.15}*),Frame->True,FrameLabel->{x,y}]

(* this may not be as bad as it looks, since the biggest effect is to add a constant to the potential, the electric field may still come out reasonably well: *)

(* h=0.01;
Ex[x_,y_]:=-(solI[x+h,y]-solI[x-h,y])/(2h) *)

(* Go for a highly accurate solution by demanding options in the finite-element approach. This will cost some CPU cycles and will gobble memory, i.e., do not expect instant return on the plots below! *)

solI=NDSolveValue[{Poisson,(Subscript[Γ, D])},V,{x,y}∈Ω,Method->{"PDEDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->0.005},"IntegrationOrder"->3}}]

(* calculate the x-component of the electric field: *)

Clear[Ex];
Ex[x_,y_]:=-D[solI[x1,y],x1]/.x1->x

Ex[0,2]

Ex0[x_,y_]:=-D[solI0[x1,y],x1]/.x1->x

PL1=Plot[Ex[x,0],{x,0,5},PlotStyle->Red];
PL2=Plot[Ex0[x,0],{x,0,5},PlotStyle->Blue];
Show[PL2,PL1,PlotRange->All]
(* the finite-element solution done to high precision yields an acceptable derived field, the original one does not. The derivative computation presented a challenge, since the potential looked OK on the contour plot. *)

Ediff[x_]:=Ex[x,0]-(Es/.s->x)

Ediff[0.1]

xT=Table[0.1j,{j,1,50}];

ListPlot[Table[{xT[[j]],Ediff[xT[[j]]]},{j,1,Length[xT]}],PlotStyle->Green,PlotRange->All]

RelError[x_]:=Abs[Ediff[x]/(Es/.s->x)]

ListPlot[100RelError[xT],PlotStyle->Green,PlotRange->All,PlotLegends->Placed[{"Relative Absolute Error [%]"},Above]]

(* this does inspire some confidence in the finite-element solution. Now solve a more interesting problem. *)

fϕ=Cos[ArcTan[x,y]]

δ=10^-9;

fϕ=x/Sqrt[x^2+y^2+δ^2];
(* without this regulator NDSolveValue complains! *)

Poisson2=D[V[x,y],{x,2}]+D[V[x,y],{y,2}]==-4Pi a^2/2 Exp[-(x^2+y^2)/a^2]fϕ

solI2=NDSolveValue[{Poisson2,(Subscript[Γ, D])},V,{x,y}∈Ω,Method->{"PDEDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->0.005},"IntegrationOrder"->4}}]

ContourPlot[solI2[x,y],{x,-x0,x0},{y,-y0,y0},ContourLabels->All,Contours->{-0.03,-0.02,-0.01,-0.005,0,0.005,0.01,0.02,0.03},Frame->True,FrameLabel->{x,y}]

Plot[solI2[x,0],{x,-x0,x0}]

(* Challenge question:

Setting the potential to zero on the disk boundary may not be the best thing to do. Thus, one should attempt an exact solution for comparison purposes!

Will the equation separate in cylindrical coordinates? *)

LaplaceCyl[g_]:=g[ϕ]D[s D[f[s],s],s]/s+f[s]D[g[ϕ],{ϕ,2}]/s^2

LaplaceCyl[Cos]

(* apparently, yes! We assume V[s,ϕ]=f[s]Cos[ϕ], divide off Cos[ϕ], and obtain an ODE for f[s]: *)

PoissonC2=D[s D[Vc2[s],s],s]/s-Vc2[s]/s^2==-4Pi a^2/2 Exp[-s^2/a^2]

solC2=DSolve[PoissonC2,Vc2[s],s]

Limit[Vc2[s]/.solC2,s->0]

Vs2=(Vc2[s]/.solC2/.{C[2]->-I C[1]}) [[1]]//FullSimplify

(* a linear rise? better, not! *)

Limit[Vs2,s->Infinity]

Vs2=Vs2/.C[1]->Pi^(3/2)/16

PL3=Plot[Vs2,{s,0,6},PlotRange->All];

PL4=Plot[solI2[x,0],{x,0,x0},PlotStyle->Red];
Show[PL3,PL4]

(* the correct procedure would be to estimate the dipole moment of the two-dimensional density Px=<s Cos[ϕ]>, and to implement a boundary condition on the basis of the correct asymptotic behavior of the potential.

Note how the numerical solution is corrupted down to values well within x0/2. *)

NetCharge=Integrate[4Pi a^2/2 Exp[-s^2/a^2]Cos[ϕ]s,{s,0,Infinity},{ϕ,0,2Pi}]

HalfCharge=Integrate[4Pi a^2/2 Exp[-s^2/a^2]Cos[ϕ]s,{s,0,Infinity},{ϕ,-Pi/2,Pi/2}]

Px=Integrate[4Pi a^2/2 Exp[-s^2/a^2]Cos[ϕ]s s Cos[ϕ],{s,0,Infinity},{ϕ,0,2Pi}]

(* It turns out that the asymptotic behaviour of the potential is ~1/s: *)

Series[Vs2,{s,Infinity,5}]//Expand

PL5=Plot[Pi^(3/2)/(128s),{s,0.6,6},PlotStyle->{Green,Dashed},PlotRange->All];

Show[PL3,PL4,PL5]

(* Note that we are comparing x-dependent solutions to s-dependence assuming ϕ=0.

The unexpected strong asymptotic behavior must be associated with the Log[s]-dependence of the potential for an infinite line of charge. *)

(* Mini-project:
Explore the option of Cos[m ϕ] dependence in the planar density for m=2, which should result in a quadrupolar planar structure. Follow the steps to derive the correct equation for f[s], and compare to the numerical solution.
 *)

(* Let us modify the boundary condition to obtain the proper numerical solution: *)

(Subscript[Γ, D2])=DirichletCondition[V[x,y]==Pi^(3/2)/(128 10a)fϕ,x^2+y^2==(10a)^2];

solI2c=NDSolveValue[{Poisson2,(Subscript[Γ, D2])},V,{x,y}∈Ω,Method->{"PDEDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->0.005},"IntegrationOrder"->4}}]

ContourPlot[solI2c[x,y],{x,-x0,x0},{y,-y0,y0},ContourLabels->All,Contours->{-0.03,-0.02,-0.01,-0.005,0,0.005,0.01,0.02,0.03},Frame->True,FrameLabel->{x,y}]

(* this looks different from the one with wrong BC! 

The funny shapes in the contours are outside the computed region, i.e., extrapolation was used to generate the data *)

PL4c=Plot[solI2c[x,0],{x,0,x0},PlotStyle->Red];

Show[PL3,PL4c,PL5]

Clear[Ex2];
Ex2[x_,y_]:=-D[solI2c[x1,y],x1]/.x1->x

PL1=Plot[Ex2[x,0],{x,δ,5},PlotStyle->Red,PlotRange->All];

Es2=-D[Vs2,s]

PL2=Plot[Es2,{s,0,5},PlotStyle->{Blue,Dashed},PlotRange->All];

Show[PL1,PL2,PlotRange->All]