function createPhaseDiagram()
    % Define axis values
    SCATTERING_LENGTH_RANGE = [79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0];
    NUM_ATOMS_LIST = [50000 60000 70000 80000];

    M = zeros(12, 4);  % Phase diagram matrix

    fig = figure('Name', 'Manual Phase Diagram Input', ...
                 'NumberTitle', 'off', ...
                 'KeyPressFcn', @keyCallback);

    hImg = imagesc(M);
    set(gca, 'YDir', 'normal');
    colormap(parula);
    colorbar;
    axis equal tight;
    xticks(1:4);
    yticks(1:12);
    xticklabels(string(NUM_ATOMS_LIST));
    yticklabels(string(SCATTERING_LENGTH_RANGE));
    xlabel('Number of Atoms');
    ylabel('Scattering Length (a0)');
    grid on;
    title('Click a cell to label it. Press "q" to quit.');

    set(fig, 'WindowButtonDownFcn', @clickCallback);
    uiwait(fig);  % Wait until user quits

    [filename, pathname] = uiputfile('phase_diagram_matrix.mat', 'Save Matrix As');
    if ischar(filename)
        save(fullfile(pathname, filename), 'M', 'SCATTERING_LENGTH_RANGE', 'NUM_ATOMS_LIST');
        disp(['Matrix saved to ' fullfile(pathname, filename)]);
    end

    function clickCallback(~, ~)
        C = get(gca, 'CurrentPoint');
        col = round(C(1,1));
        row = round(C(1,2));
        if row >= 1 && row <= 12 && col >= 1 && col <= 4
            prompt = sprintf('Enter label for scattering length %.1f and N = %d:', ...
                             SCATTERING_LENGTH_RANGE(row), NUM_ATOMS_LIST(col));
            answer = inputdlg(prompt, 'Input Label', 1);
            if ~isempty(answer)
                val = str2double(answer{1});
                if ~isnan(val)
                    M(row, col) = val;
                    refreshDisplay();
                else
                    warndlg('Invalid input. Please enter a numeric value.', 'Invalid Input');
                end
            end
        end
    end

    function keyCallback(~, event)
        if strcmp(event.Key, 'q')
            uiresume(fig);
            close(fig);
        end
    end

    function refreshDisplay()
        set(hImg, 'CData', M);
        drawnow;
    end
end