Построение двумерного уравнения диффузии во времени

Question

Итак, в основном, я создал график на python, который моделирует две взаимодействующие популяции на острове и показывает использование уравнения диффузии для моделирования движения и изменения популяции в одном измерении, и я пытаюсь изменить еготак что вместо этого у меня есть двухмерный график, который учитывает положение х и у видов на острове.

Извиняюсь за то, что код довольно длинный, хотя большая его часть просто определяет параметры - я не слишком уверен, как создать двухмерный график с плотностью населения в каждой позиции, это мойпопытка ниже, но сюжет не является правильным.Я думаю, что мое численное решение должно быть правильным, я предполагаю, что это мои начальные / граничные условия и как я настроил свой массив для распределения вероятностей, что является неправильным, но я не уверен, как его изменить.Например, одна из проблем, с которыми я сталкиваюсь, заключается в том, что она имеет позицию Y на оси x - как я могу это исправить?

r1=0.1 # growth rate B
r2=0.3 # growth rate M
KB = 20 # carrying capacity B
KM = 6000 # carrying capacity M
x=np.arange(0,70) # distance from the coast (100m)
y=np.arange(0,90) # distance from the coast
dx=1 # change in x position
dy=1 # change in y position
dt=1 # time step
m=71 # number of x steps
o=91 # number of y steps
n=21 # time
alpha =  0.00002 # predation rate
beta = 0 # growth rate from predation (assumed to be zero)

B=np.zeros(shape=(m,o,n)) # species B 
M=np.zeros(shape=(m,o,n)) # species M
D=0.35 # diffusivity of B
D2=0.05 # diffusivity of M
Alpha=(D*dt)/(dx*dx)
Beta=(D*dt)/(dy*dy) 
Alpha2=(D2*dt)/(dx*dx)
Beta2 = (D2*dt)/(dy*dy)

M[0,:,0]=0 #initial conditions
M[m-1,:,0]=0
M[:,0,0]=0
M[:,m-1,0]=0
B[1:26,:,0]=0.96
B[26:44,:,0]=4.6
B[44:m-1,:,0]=0.96
B[:,1:26,0]=0.96
B[:,26:44,0]=4.6
B[:,44:m-1,0]=0.96
M[1:26,:,0]=2500
M[26:44,:,0]=1000
M[44:m-1,:,0]=2500
M[:,1:26,0]=2500
M[:,26:44,0]=1000
M[:,44:m-1,0]=2500
for k in range(0,n-1): # loop for time
    B[0,:,k+1]=B[1,:,k+1] # boundary conditions
    B[m-1,:,k+1]=B[m-2,:,k+1]
    M[0,:,k+1]=0
    M[m-1,:,k+1]=0
    B[:,0,k+1]=B[:,1,k+1]
    B[:,m-1,k+1]=B[:,m-2,k+1]
    M[:,0,k+1]=0
    M[:,m-1,k+1]=0
    for i in range(1,m-1): # loop for x
        for j in range(1,o-1): # loop for y
            Bdxx=(-2*Alpha)*B[i,j,k]+Alpha*B[i+1,j,k]+Alpha*B[i-1,j,k]
            Bdyy=(-2*Beta)*B[i,j,k]+Beta*B[i,j+1,k]+Beta*B[i,j-1,k]
            B[i,j,k+1]=B[i,j,k]+Bdxx+Bdyy+r1*B[i,j,k]*(1-B[i,j,k]/KB)-alpha*M[i,j,k]*B[i,j,k] # numerical solution for B 
            Mdxx=(-2*Alpha2)*M[i,j,k]+Alpha2*M[i+1,j,k]+Alpha2*M[i-1,j,k]
            Mdyy=(-2*Beta2)*M[i,j,k]+Beta2*M[i,j+1,k]+Beta2*M[i,j-1,k]
            M[i,j,k+1]=M[i,j,k]+Mdxx+Mdyy+r2*M[i,j,k]*(1-M[i,j,k]/KM)+beta*M[i,j,k]*B[i,j,k] # numerical solution for M


plt.pcolormesh(B[:,:,n-1]) 
plt.colorbar()
plt.xlabel('X Position')
plt.ylabel('Y Position')

Любая помощь действительно приветствуется, спасибо!

редактировать: я понял, что часть моей проблемы заключается в том, что начальные условия для координат у неправильные