i am doing finite difference scheme to solve Fisher kpp equation :$u_{t}=u_{xx}+u_{yy}+u-u^2 and with initial condition u(x,y,0)=1 for (x,y)=(0,0)and 0 other wise and u=0 at boundary ,also plot showing the positions where u(x, y, t) = 0.5 for each time step .

here the Python code :

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps):
    u = np.zeros((Nx, Ny))
    u[0, 0] = 1  # Set initial condition at the center

    Dxx, Dyy = 1/dx**2, 1/dy**2
    i = slice(1, Nx-1)
    iw = slice(0, Nx-2)
    ie = slice(2, Nx)
    j = slice(1, Ny-1)
    js = slice(0, Ny-2)
    jn = slice(2, Ny)

    for _ in range(steps):
        dux = Dxx * (u[ie, j] - 2*u[i, j] + u[iw, j])
        duy = Dyy * (u[i, jn] - 2*u[i, j] + u[i, js])
        reaction = r * u[i, j] * (1 - u[i, j])
        u[i, j] += dt * (D * (dux + duy) + reaction)

        # Set boundary conditions
        u[0, :] = 0
        u[-1, :] = 0
        u[:, 0] = 0
        u[:, -1] = 0

    return u

# Parameters
Nx = 100        # Number of mesh points in the x-direction
Ny = 100        # Number of mesh points in the y-direction
D = 1.0         # Diffusion coefficient
r = 1.0         # Reaction rate
dx = 1 / Nx
dy = 1 / Ny
dt = (dx**2 * dy**2) / (2 * D * (dx**2 + dy**2))

# Time steps for t = 5, t = 10, t = 20
steps_t5 = int(5 / dt)
steps_t10 = int(10 / dt)
steps_t20 = int(20 / dt)

# Solve the Fisher-KPP equation for t = 5, t = 10, and t = 20
solution_t5 = solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps_t5)
solution_t10 = solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps_t10)
solution_t20 = solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps_t20)

# Create meshgrid for plotting
x = np.linspace(0, 1, Nx)
y = np.linspace(0, 1, Ny)
X, Y = np.meshgrid(x, y)
# Collecting positions where u(x, y, t) = 0.5 for each time step
T = 10
desired_positions = []

for t in range(1, T + 1):
    steps_T = int(t / dt)
    solution_T = solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps_T)
    
    # Finding (x, y) such that u(x, y, t) = 0.5 using np.where
    half_max = np.max(solution_T) / 2
    desired_idx_x, desired_idx_y = np.where(solution_T == half_max)
    desired_x = desired_idx_x * dx
    
    desired_positions.append((t, desired_x))

# Plotting positions x as a function of time step
fig, ax = plt.subplots()
for t, x_positions in desired_positions:
    ax.plot([t] * len(x_positions), x_positions, 'ro')

ax.set_xlabel('Time Step')
ax.set_ylabel('Position x')
ax.set_title('Positions where u(x, y, t) = 0.5')
plt.show()

I’m expecting that the graph slope of (the graph of positions in function of time )to be approximatively equal or large then2 but my code doesnt run ,where is the issue?