import numpy
import pylab
from mpl_toolkits.mplot3d import Axes3D

slices = 100.

bound = 10.

width = 2.*bound/slices

enot = 8.85*10**-12

rho = (1.5*10**-7)/(8*bound**3)

x = numpy.arange(0,2*bound, width)
y = numpy.arange(0,2*bound, width)
z = numpy.arange(0,2*bound, width)

xx = numpy.zeros(x.size)
yy = numpy.zeros(y.size)
zz = numpy.zeros(z.size)

for i in range(0,x.size):
    xx[i] = -bound + i*width
    yy[i] = -bound + i*width
    zz[i] = -bound + i*width

v = numpy.zeros([100,100,100])

v[0,:,:] = xx
v[:,0,:] = yy
v[:,:,0] = zz

change = 1.
while change < 10:
    for i in range(1, xx.size-1):
            for j in range(1, xx.size-1):
                for k in range(1, xx.size-1):
                    v[i,j,k] = (1./6)*(v[i+1,j,k]+v[i-1,j,k] + v[i,j+1,k]+v[i,j-1,k] + v[i,j,k+1]+v[i,j,k-1]) + (width**2/(6*enot))*rho
    print change
    change = change + 1

fig = pylab.figure()
ax = Axes3D(fig)
ax.scatter(v[0,:,:], v[:,0,:], v[:,:,0])
pylab.show()