To visualize the solution of the foolowing pde :$u_t=u_xx+u-yy+f(u)$at different times, such as

t=1 and t=3, i use finite diffrence scheme and the following Python code :

this the dode i use :

import numpy as np
import matplotlib.pyplot as plt

# Parameters
D = 0.1  # Diffusion coefficient
L = 1.0  # Length of domain
T = 3.0  # Total time
Nx = 100  # Number of grid points in x-direction
Ny = 100  # Number of grid points in y-direction
Nt = 3000  # Number of time steps
dx = L / (Nx - 1)  # Grid spacing in x-direction
dy = L / (Ny - 1)  # Grid spacing in y-direction
dt = T / Nt  # Time step size

# Initial condition
def initial_condition(x, y):
    return np.exp(-((x - 0.5)**2 + (y - 0.5)**2) / 0.1)

# Reaction term
def reaction_term(u):
    return u * (1 - u)

# Set boundary condition (u = 0 at boundary)
def apply_boundary_condition(u):
    u[0, :] = 0  # Bottom boundary
    u[-1, :] = 0  # Top boundary
    u[:, 0] = 0  # Left boundary
    u[:, -1] = 0  # Right boundary

# Create grid
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# Initialize solution array
u = np.zeros((Nx, Ny))

# Set initial condition
u = initial_condition(X, Y)

# Time points to visualize
times_to_visualize = [1, 3]  # Time points at which to visualize solution

# Iterate over time
for n in range(Nt + 1):
    # Compute spatial derivatives using central differences
    u_xx = (np.roll(u, -1, axis=0) - 2*u + np.roll(u, 1, axis=0)) / dx**2
    u_yy = (np.roll(u, -1, axis=1) - 2*u + np.roll(u, 1, axis=1)) / dy**2
    
    # Compute reaction term
    f_u = reaction_term(u)
    
    # Update solution using forward Euler method
    u += dt * (D * (u_xx + u_yy) + f_u)
    
    # Apply boundary condition
    apply_boundary_condition(u)
    
    # Check if current time is in times_to_visualize
    if n * dt in times_to_visualize:
        # Plot solution
        plt.figure()
        plt.contourf(X, Y, u, cmap='coolwarm')
        plt.colorbar(label='u')
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title('2D Reaction-Diffusion Equation at t = {:.2f}'.format(n * dt))
        plt.show()

and mu aim is to have something like this :
enter image description here
but my code doesnt give me such result but a strange result here what i getenter image description here
where im mistaken ? and how to correct my code ( im not usully person who do a lot of numerical analysis so i’m sorry if it seems that the question is not advanced )