The brine tank is a historical and excellent way to find solutions to the Laplace equation! This was used to design things like cyclotron magnets and electrostatic lenses for charged particles before computers were fast enough.
Image from Tanabe and Yamada 1958; Science reports of the Research Institute, Tohoku University, Series A, 10, 133-174
The simplest way to solve it numerically is to do "Jacobi relaxation". Read more about it here for example. In cartesian coordinates, you make a numerical grid and set the points that correspond to the surface of the conductors to their potential, and set the potential to the boundary to some intermediate value. Below I've chosen -1, +1 and 0. The boundary needs to be far enough away that it doesn't affect your region of interest much. This is an approximate technique.
Then you iterate. For each pass, you make a new potential that is the average of the +/-x and +/-y nearest neighbors, four in total. You update only the points that are not on the conductor or boundary. After thousands of iterations the potential approaches a solution to the Laplace equation.
If you want a better solution, make the spacing smaller and the boundary farther away. If it's too slow (frequently it is in 3D) then you can use
Multigrids; solve on a coarse (fast) grid, then interpolate to a fine grid and iterate a little longer
Successive Overrelaxation; "overshoot" the averaging by changing the values by an amount larger than 100% of what you would have with simple averaging. Try values like 110% to 150%, it can go unstable if they are too large.
If you have any questions or need clarification, leave a comment. I used this technique in this question, where I saved intermediate snapshots and turned them into a GIF.
Here is an example in Python:
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import convolve
d = 0.25
x = np.arange(0, 35+d, d)
y = np.arange(0, 48+d, d)
X, Y = np.meshgrid(x, y)
phi = np.zeros_like(X)
isneg = np.zeros_like(X, dtype=bool)
ispos = np.zeros_like(X, dtype=bool)
bound = np.zeros_like(X, dtype=bool)
isneg[(X>=9)*(X<=19)*(Y>=10)*(Y<=12)] = True
ispos[(X>=10)*(X<=20)*(Y>=28)*(Y<=30)] = True
ispos[(X>=18)*(X<=20)*(Y>=20)*(Y<=28)] = True
bound[:,0] = bound[:,-1] = bound[0,:] = bound[-1,:] = True
phi[isneg] = -1.0
phi[ispos] = +1.0
phi[bound] = 0.0
updateme = np.ones_like(X, dtype=bool)
updateme[isneg] = updateme[ispos] = updateme[bound] = False
kernel = 0.25*np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]], dtype=float)
phi0 = phi.copy() # keep the original handy
keepers = 0, 100, 500, 1000, 2000, 9999
keep = 
for i in range(10000):
if i in keepers:
phi2 = convolve(phi, kernel)
phi[updateme] = phi2[updateme]
for i, (thing, n) in enumerate(zip(keep, keepers)):
plt.subplot(2, 3, i+1)