from pylab import *

m = 2.0		# Mass of the particle
q = 5.0		# Charge
t = 0.0
dt = 0.001
N = 5000	# number of integrating steps

E = array([0.0, 0.0, 0.0])
B = array([0.0, 0.0, 5.0])

pos = zeros([N,3])   # array of N positions
pos[0] = [5,0,0]     # inital position. x, y and z
vel = zeros([N,3])
vel[0] = [20,0,0]    # initial velocity. vx, vy and vz

for k in range(N-1):
	F = q * (E + cross(vel[k],B)) 		# F = q * [E (v x B)]
	vel[k+1] = vel[k] + F/m * dt		# a = F/m; dv = a.dt
	pos[k+1] = pos[k] + vel[k] * dt		# dx = v.dt
	t = t + dt

from mpl_toolkits.mplot3d import Axes3D
ax = Axes3D(figure())
ax.plot(pos[:,0], pos[:,1], pos[:,2])		# 3d plot of x, y and z
ax.set_zlabel('Z axis')
show()