% Parameters
L = 10;                  % Length of domain along x
T = 500;                  % Total time duration
nx = 100;                % Number of grid points along x
nt = 5000;               % Number of time steps
dx = L/nx;               % Grid size along x
dt = T/nt;               % Time step size
b=0.01;
% Reaction terms parameters
r = 0.1;
K = 10;
a = 0.03;

% Diffusion coefficients
D = 0.02;
C = 0.01;

% Initialize the grid values
Y = zeros(nx+1, nt+1);
X = zeros(nx+1, nt+1);

% Initial conditions
Y(:,1) = 2
X(:,1) = 1
% Time stepping loop
for t = 1:nt               
    % Interior points
    for i = 2:nx
        d2Y = (Y(i + 1, t) - 2 * Y(i, t) + Y(i - 1, t)) / dx^2;
        d2X = (X(i + 1, t) - 2 * X(i, t) + X(i - 1, t)) / dx^2;
        
        Y(i, t + 1) = Y(i, t) + dt * (r * Y(i, t) * X(i,t) * (1 - Y(i, t) / K) + D * d2Y);
        X(i, t + 1) = X(i, t) + dt * (a * Y(i, t) - b*X(i,t).^2 + C * d2X);
    end
    
    % Boundary conditions
    Y([1, nx + 1], t + 1) = Y([1, nx + 1], t);
    X([1, nx + 1], t + 1) = X([1, nx + 1], t);
end

% Plot results
t = (0:nt)*dt;
x = (0:nx)*dx;

figure;
subplot(1,2,1);
mesh(t, x, Y);
xlabel('Time');
ylabel('Position');
zlabel('Population');
title('opulation vs Time & Position');

subplot(1,2,2);
mesh(t, x, X);
xlabel('Time');
ylabel('Position');
zlabel('Information Density');
title('Information Density vs Time & Position');

figure;
plot(X(:,end), Y(:,end));
xlabel('iiPopulation Density');
ylabel('ppInformation Density');
title('Information Density vs Population Density');
% Sum over the position for each time step
total_population = sum(Y, 1);

% Plot total population against time
figure;
plot(t, total_population);
xlabel('Time');
ylabel('Total Population');
title('Total Population vs Time');
what parameters to change to slow down the growth of the total population