Generate figure 2 of paper Role of dissolved oxygen on the plankton dynamics in Spatio-temporal domain

Question

I am trying to generate figure 2 of the above paper.

I generated the following code with the help of ChatGPT.

% Parameters
r = 0.445532;
a = 3;
eta = 0.2;
delta1 = 0.1;
beta = 0.2;
beta1 = 0.4;
D0 = 3;
alpha1 = 0.18;
alpha2 = 0.2;
delta2 = 0.3;
D_a = 0.06;
D_b = 0.04;
D_c = 0.02;

% Spatial and time domain
L = 10;
N = 50; % Number of spatial points (adjust as needed)
dx = L / (N - 1);
dt = 0.1; % Time step (adjust as needed)
T = 300;
M = round(T / dt);

% Initialize variables
P = ones(N, N);
D = ones(N, N);
Z = ones(N, N);

% Boundary conditions
P(:, 1) = P(:, 2);
P(:, N) = P(:, N-1);
P(1, :) = P(2, :);
P(N, :) = P(N-1, :);

D(:, 1) = D(:, 2);
D(:, N) = D(:, N-1);
D(1, :) = D(2, :);
D(N, :) = D(N-1, :);

Z(:, 1) = Z(:, 2);
Z(:, N) = Z(:, N-1);
Z(1, :) = Z(2, :);
Z(N, :) = Z(N-1, :);

% Finite difference method
for m = 1:M
    for i = 2:N-1
        for j = 2:N-1
            laplacian_P = (P(i+1,j) + P(i-1,j) + P(i,j+1) + P(i,j-1) - 4 * P(i,j)) / dx^2;
            laplacian_D = (D(i+1,j) + D(i-1,j) + D(i,j+1) + D(i,j-1) - 4 * D(i,j)) / dx^2;
            laplacian_Z = (Z(i+1,j) + Z(i-1,j) + Z(i,j+1) + Z(i,j-1) - 4 * Z(i,j)) / dx^2;
            
            dP_dt = (r * P(i,j) / (D0 - D(i,j)) - delta1 * P(i,j) - (beta * P(i,j) * Z(i,j)) / (P(i,j) + a) + D_a * laplacian_P) * dt;
            dD_dt = (eta * (D0 - D(i,j)) - alpha1 * P(i,j) - alpha2 * Z(i,j) + D_b * laplacian_D) * dt;
            dZ_dt = (beta1 * P(i,j) * Z(i,j) / (P(i,j) + a) - delta2 * Z(i,j) + D_c * laplacian_Z) * dt;
            
            P(i,j) = P(i,j) + dP_dt;
            D(i,j) = D(i,j) + dD_dt;
            Z(i,j) = Z(i,j) + dZ_dt;
        end
    end
end

% Plot results
[X, Y] = meshgrid(1:N, 1:N);

figure;
subplot(1, 3, 1);
contourf(X, Y, P, 20, 'LineStyle', 'none');
title('P');
colorbar;

subplot(1, 3, 2);
contourf(X, Y, D, 20, 'LineStyle', 'none');
title('D');
colorbar;

subplot(1, 3, 3);
contourf(X, Y, Z, 20, 'LineStyle', 'none');
title('Z');
colorbar;

sgtitle('Plankton Dynamics at t = 300');

The following is the result.

I understand that the boundary conditions and initial conditions I use seem not correct. Any help is much appreciated.