python - fast method with numpy for 2D Heat equation -

i'm looking method solve 2d heat equation python. have implemented finite difference method slow motion (to make 100,000 simulations takes 30 minutes). idea create code in end can write,

for t in time:       deltau=f(u)     u=u+deltau*deltat     save(u) 

how can that?

in first form of code, used 2d method of finite difference, grill 5000x250 (x, y). decrease speed of computing , idea find

deltau = f(u)

where u heat function. implementation used source http://www.timteatro.net/2010/10/29/performance-python-solving-the-2d-diffusion-equation-with-numpy/ 2d case, run time more expensive necessity. there methods this?

maybe must work matrix

a=1/dx^2 (2 -1  0  0 ... 0           -1  2 -1  0 ... 0          0  -1  2 -1 ... 0          .            .          .            .          .            .          0 ...        -1 2) 

but how make in 2d problem? how inserting boundary conditions in a? code finite difference used:

lx=5000   # physical length x vector in micron ly=250   # physical length y vector in micron nx = 100     # number of point of mesh along x direction ny = 50     # number of point of mesh along y direction = 0.001 # diffusion coefficent dx = 1/nx dy = 1/ny dt = (dx**2*dy**2)/(2*a*(dx**2 + dy**2)) # 0.04 x = linspace(0.1,lx, nx)[np.newaxis] # vector create mesh y = linspace(0.1,ly, ny)[np.newaxis] # vector create mesh i=sqrt(x*y.t) #initial data heat equation u=np.ones(([nx,ny])) # u matrix referred heat function steps=100000 m in range (0,steps):         du=np.zeros(([nx,ny]))          in range (1,nx-1):              j in range(1,ny-1):                 dux = ( u[i+1,j] - 2*u[i,j] + u[i-1, j] ) / dx**2                duy = ( u[i,j+1] - 2*u[i,j] + u[i, j-1] ) / dy**2                            du[i,j] = dt*a*(dux+duy)         # boundary conditions     t1=(u[:,0]+u[:,1])/2     u[:,0]=t1     u[:,1]=t1     t2=(u[0,:]+u[1,:])/2     u[0,:]=t2     u[1,:]=t2     t3=(u[-1,:]+u[-2,:])/2     u[-1,:]=t3     u[-2,:]=t3     u[:,-1]=1       filename1='data_{:08d}.txt'      if m%100==0:         np.savetxt(filename1.format(m),u,delimiter='\t' ) 

for elaborate 100000 steps run time 30 minutes. optimize code (with idea presented in initial lines) have run time 5/10 minutes or minus. how can it?

have considered paralellizing code or using gpu acceleration.

it if ran code python profiler (cprofile) can figure out bottleneck in runtime is. i'm assuming it's in solving matrix equation can sped methods listed above.