Simulate / plot electrostatic field I have something like this: 

For given voltage (9 V red wire, 0 V blue wire) I want to simulate electrostatic field. 
I already know, for example that in the point (10.5;0) there are 10 V and etc.
I want to simulate and plot electrostatic field for future comparison and for my report.
Is there any software for doing that? I like to work with Maple, so if there is any toolboxes, I would be glad to hear about them.
For slightly different position plot looks awful now: 

 A: As I said earlier in my comment, Comsol is pretty useful for these sort of things. I took your problem, and this is the plot that I got. I hope it's useful. 
 
A: 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:
        keep.append(phi.copy())
        print i
    phi2 = convolve(phi, kernel)
    phi[updateme] = phi2[updateme]

plt.figure()
for i, (thing, n) in enumerate(zip(keep, keepers)):
    plt.subplot(2, 3, i+1)
    plt.imshow(thing, origin='lower')
    plt.title("n="+str(n))
plt.show()

