Generating Turing Patterns Using Finite Difference on Reaction Diffusion Equations

Question

For a school project, I am implementing in Python a finite difference method to numerically solve the following system of reaction diffusion PDE:

I have so far implemented the following:

import numpy as np
import matplotlib.pyplot as plt

D_u = 1.0
D_v = .03
k1 = 5.0
k2 = 11.0
dx = .01
dy = .01
dt = .0001 

# Initialize domain
x = np.arange(0, 3, dx)
y = np.arange(0, 3, dy)

plate_length = len(x)
max_iter_time = 10000

# initialize solution matrix
u = np.zeros((max_iter_time, plate_length, plate_length)) #u(t,x,y) 
v = np.zeros((max_iter_time, plate_length, plate_length)) #v(t,x,y)

# Initial condition : (u,v) = (u_s, v_s) + r
random_perturbation_u = np.random.normal(loc = 0.01, size=(plate_length, plate_length))
random_perturbation_v = np.random.normal(loc = 0.01, size=(plate_length, plate_length))

u[0,:,:] = 1 + (k2**2)/25 + random_perturbation_u
v[0,:,:] = k2/5 + random_perturbation_v

def periodic_bc(matrix):
    matrix[:, 0, :] = matrix[:, -2, :] 
    matrix[:, -1, :] = matrix[:, 1, :]  
    matrix[:, :, 0] = matrix[:, :, -2]  
    matrix[:, :, -1] = matrix[:, :, 1]  

    return matrix

#implement finite difference
def calculate_u_and_v(u,v):
    for k in range(0, max_iter_time-1,1):
        for j in range(1, plate_length-1):
            for i in range(1, plate_length-1): 
                u[k + 1, i, j] = u[k][i][j] + dt * ((D_u * ((u[k][i+1][j] - 2*u[k][i][j]+ u[k][i-1][j])/dx**2 + (u[k][i][j+1] - 2*u[k][i][j] + u[k][i][j-1])/dy**2))
                + k1 * (v[k][i][j] - (u[k][i][j] * v[k][i][j])/(1 + v[k][i][j]**2)))

                v[k + 1, i, j] = v[k][i][j] + dt * ((D_v * ((v[k][i+1][j] - 2*v[k][i][j]+ v[k][i-1][j])/dx**2 + (v[k][i][j+1] - 2*v[k][i][j] + v[k][i][j-1])/dy**2))
                + k2 - v[k][i][j] - (4 * u[k][j][i] * v[k][j][i])/(1 + v[k][i][j]**2))

        u = periodic_bc(u)
        v = periodic_bc(v)
        
    return u,v

u,v = calculate_u_and_v(u,v)

# plot the final image
plt.figure()
plt.pcolormesh(x, y, v[max_iter_time-1], cmap=plt.cm.jet, vmin=np.min(v), vmax=np.max(v))
plt.title(f"t = {max_iter_time}")
plt.ylabel(f"$k_1=${k1}")
plt.xlabel(f"$D_v=${D_v}")
plt.colorbar(ticks = [])

plt.show()

From this I get results like this: my results

But I am trying to replicate/reproduce patterns like this: desired result

Does anyone have any idea where my implementation goes wrong? I am still relatively new to this type of thing. Thank you very much.