<<"VectorAnalysis`"

CoordinateSystem

Coordinates[]

SetCoordinates[Cartesian[x,y,z]]

Grad[x^2 +y^2]

ϕ[x,y,z]:=(1/(Sqrt[(x^2)+(y^2)+(z^2)]))

e:=-Grad[ϕ[x,y,z]]

e

Div[e]

Curl[e]

ϕ1[Q_,x0_,y0_,z0_]:=(Q/(4 π ϵ0 (Sqrt[(((x-x0))^2)+(((y-y0))^2)+(((z-z0))^2)])))

ϕ2 = ϕ1[Q,0,0,a] + ϕ1[-Q,0,0,-a]

ContourPlot[Evaluate[ϕ2/.x->0/.Q->1/.ϵ0->1/.a->1/2],{y,-1,1},{z,-1,1},Contours->40,PlotPoints->50,ColorFunction->(Hue[1/2-#/2]&)]

DensityPlot[Evaluate[ϕ2/.x->0/.Q->1/.ϵ0->1/.a->1/2],{y,-2,2},{z,-2,2},PlotPoints->200,Mesh->False,ColorFunction->(Hue[1-#/2]&)]

ϕ4 = ϕ1[Q,0,-a,a] + ϕ1[-Q,0,-a,-a]+ϕ1[Q,0,a,-a]+ϕ1[-Q,0,a,a]

ContourPlot[Evaluate[ϕ4/.x->0/.Q->1/.ϵ0->1/.a->1/2],{y,-1,1},{z,-1,1},Contours->40,PlotPoints->50,ColorFunction->(Hue[1-#/2]&)]

ϕ5[Q_,x_,y_,z_,x0_,y0_,z0_]:=(Q/(4 π ϵ0 (Sqrt[(((x-x0))^2)+(((y-y0))^2)+(((z-z0))^2)])))

e5[Q_,x_,y_,z_,x0_,y0_,z0_]:=-Grad[ϕ5[Q,x,y,z,x0,y0,z0]]

e=e5[Q,x,y,z,x0,y0,z0]