function GrayScottModel
    % Parameters
    da = 0.2;    % Diffusion rate of u
    db = 0.1;    % Diffusion rate of v
    f = 0.035;   % Feed rate
    g = 0.06;    % Kill rate
    k = 1;       % Reaction rate

    % Grid size and initialization
    grid_size = 256;
    u = ones(grid_size, grid_size);
    v = zeros(grid_size, grid_size);

    % Seed small area with v concentration
    mid = grid_size / 2;
    range = mid-10:mid+10;
    v(range, range) = 0.6;

    % Create figure and set up mouse click callback
    hFig = figure('Name', 'Gray-Scott Model', 'NumberTitle', 'off');
    hAx = axes('Parent', hFig);
    set(hFig, 'WindowButtonDownFcn', @mouseClick);

    % Time-step settings
    dt = 1;
    n_steps = 5000;

    for t = 1:n_steps
        % Compute Laplacians
        Lu = del2(u);
        Lv = del2(v);

        % Reaction-diffusion equations
        u_new = u + dt * (da * Lu - k* u .* v.^2 + f * (1 - u));
        v_new = v + dt * (db * Lv + k * u .* v.^2 - (f + g) * v);

        u = u_new;
        v = v_new;

        % Update plot every 10 steps
        if mod(t, 10) == 0
            imagesc(hAx, v);
            axis image;
            title(sprintf('Time step: %d', t));
            colorbar;
            caxis([0 1]);  % Ensure colorbar range stays fixed from 0 to 0.5
            drawnow;
        end
    end

    function mouseClick(~, ~)
        coords = get(hAx, 'CurrentPoint');  % Get mouse coordinates
        x = round(coords(1, 1));
        y = round(coords(1, 2));
        if x > 0 && x <= grid_size && y > 0 && y <= grid_size
            v(max(1, y-2):min(grid_size, y+2), max(1, x-2):min(grid_size, x+2)) = 0.6;
        end
    end
end

function output = del2(input)
    % Discrete Laplacian in 2D
    %output = conv2(input, [0.25, 0.5, 0.25; 0.5, -3, 0.5; 0.25, 0.5, 0.25], 'same');
    output = conv2(input,[0 1 0; 1 -4 1; 0 1 0], 'same');
end
