Laplace equation

Numerical solution by an iterative  (relaxation) method on a grid in two dimensions.
Modified SolveLaplace2.nb to immerse the electrodes inside the tank, and to use an insulated metal tank boundary (current will not flow into it). 
This is achieved by demanding a vanishing normal derivative of the potential at the tank boundaries.

set up the grid and the potential on the boundaries

We solve a mathematical problem and use dimensionless variables. We have an elecrode at high potential, and another one at zero. We then choose a grid spacing h for x and y, and a number of points, from which we derive the scale (width W and depth Y)

Va=10;
Vb=-10;

h=0.1;
nx=120;
ny=60;

W=nx h;
Y=ny h;
Yh=Y/2;
Hh=Yh/4;
X1=W/4+Hh;
X2=3W/4-Hh;

We apply a simple boundary condition on the perfectly conducting outer wall, namely V=0.

We set up a table to hold the potential (Phi) and an index table (iP) which decides whether the values will be iterated or kept fixed (boundary values)

Phi=Table[0,{i,0,nx},{j,0,ny}];
iP=Table[0,{i,0,nx},{j,0,ny}];

Now set the potential at the electrodes themselves (solid cylinder - or just shell):

Do[
x=(i1-1)h;
y=(j1-1)h;
If[(x-X1)^2+(y-Yh)^2≤Hh^2,
Phi=ReplacePart[Phi,{i1,j1}->Va];
iP=ReplacePart[iP,{i1,j1}->1];
]
If[(x-X2)^2+(y-Yh)^2≤Hh^2,
Phi=ReplacePart[Phi,{i1,j1}->Vb];
iP=ReplacePart[iP,{i1,j1}->1];
]
,{i1,1,nx+1},{j1,1,ny+1}]

The complete boundary conditions are set now.
We define an iteration procedure Sweep. It steps through the array, identifies from iP whether the element should be updated (iP = 0), and performs an update by trying to enforce locally the Laplace equation in finite-difference representation.
A Jacobi iteration would keep all entries from the previous iteration, and update the array only after the sweep is completed.
A Gauss-Seidel iteration changes the Phi entries on the fly.
In the beginning the potential is zero everywhere, except on the boundaries.

The condition of vanishing normal derivative at the boundary will be implemented by re-setting the boundary values to the innter neighbor values after each sweep.

Now carry out a large number of sweeps (you should check how converged the values are after 2000 complete sweeps!):

Do[
Do[
If[iP[[i,j]]==0,
Phi[[i,j]]=(Phi[[i+1,j]]+Phi[[i-1,j]]+Phi[[i,j+1]]+Phi[[i,j-1]])/4.;
]
,{i,2,nx},{j,2,ny}]
Do[
y=(j1-1)h;
Phi=ReplacePart[Phi,{1,j1}->Phi[[2,j1]]];
Phi=ReplacePart[Phi,{nx+1,j1}->Phi[[nx,j1]]];
,{j1,1,ny+1}]
Do[
x=(i1-1)h;
Phi=ReplacePart[Phi,{i1,1}->Phi[[i1,2]]];
Phi=ReplacePart[Phi,{i1,ny+1}->Phi[[i1,ny]]];
,{i1,1,nx+1}]

,{k,1,15000}]

ListContourPlot[Phi,InterpolationOrder->3,Contours->{-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10},ContourLabels->All]

It takes over 1000 iterations to see solutions that respect the expected symmetries. Be critical of the results obtained!

So far we plotted the list Phi over the matrix indices i,j. We really need to connect {i,j} to {x_i,x_j}!

In the next step we want to compute an interpolation of the potential so that we can compute the electric field.

fT=Table[{{(i-1)h,(j-1)h},Phi[[i,j]]},{i,1,nx+1},{j,1,ny+1}];

fT[[2,2]]

fV=Interpolation[Flatten[fT,1]]

By default this used a third-order spline-type interpolation.

fV[2.,3.]

ContourPlot[fV[x,y],{x,0,W},{y,0,Y},Contours->{-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10},ContourLabels->All,AspectRatio->Automatic]

Plot3D[fV[x,y],{x,0,W},{y,0,Y},Mesh->None,ColorFunction->Hue]

d=h/2;

Define the two-dimensional electric field as the (negative) gradient of the potential:

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

Ey[x_,y_]:=-(fV[x,y+h]-fV[x,y-h])/(2h);

VectorPlot[{Ex[x,y],Ey[x,y]},{x,0,W},{y,0,Y},AspectRatio->Automatic]

StreamDensityPlot[{Ex[x,y],Ey[x,y]},{x,h,W-h},{y,h,Y-h},ColorFunction->"Rainbow",AspectRatio->Automatic]

The plot of the electric field gives us also the current density plot.
The color coding tells us where the current density is large (compare with the simple vector arrow plot for E-field)

The calculation assumes an infinite extent in the perpendicular direction (no 2nd derivatives with respect to z).

Evidently, the equipotential lines are normal to the container walls, therefore preventing charge flow onto the walls.
The potential drop along the wall is no longer linear (as it was in the different example SolveLaplace1.mw), but in the region between the electrodes along the sidewalls it is approximately so.