function plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, titleString)
    % Region labels and corresponding colors
    regionNames = ["Unmodulated", "SSD", "Stripe", "Labyrinth", "Honeycomb"];
    regionColors = [
        0.8 0.8 0.8;   % unmodulated - light gray
        0.2 0.6 1.0;   % SSD - blue
        0.2 0.8 0.2;   % Stripe - green
        1.0 0.6 0.2;   % Labyrinth - orange
        0.8 0.2 0.8    % Honeycomb - purple
    ];

    % Interpolate the discrete-valued matrix using 'nearest' for smoothing
    [X, Y] = meshgrid(NUM_ATOMS_LIST, SCATTERING_LENGTH_RANGE);
    [Xq, Yq] = meshgrid(linspace(min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST), 500), ...
                        linspace(min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE), 500));
    Mq = interp2(X, Y, M, Xq, Yq, 'nearest');

    % Plot
    figure('Color','w');  % Correctly assign figure handle
    clf
    set(gcf, 'Position', [100, 100, 950, 800]);
    set(gcf, 'Renderer', 'opengl');
    imagesc([min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST)], ...
            [min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE)], ...
            Mq);
    % Overlay smooth contour lines to show boundaries
    hold on;
    contour(Xq, Yq, Mq, 0.5:1:3.5, 'k', 'LineWidth', 1.2);  % 0.5 to 3.5 separates 0–1, 1–2, etc.
    set(gca, 'YDir', 'normal');  % Keep Y increasing upward
    colormap(regionColors);
    cb = colorbar('Ticks', 0:4, 'TickLabels', regionNames, 'FontSize', 12);
    cb.Color = 'k';  % set colorbar text color
    clim([0 4]);  % Match value range to colormap rows
    % Axis formatting
    xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
    ylabel('Scattering Length (\times a_o)', 'Interpreter', 'tex', 'FontSize', 16);
    t = title(titleString, 'Interpreter', 'tex', 'FontSize', 18);
    t.Color = 'k';  % force title color to black

    set(gca, 'FontSize', 16, ...
        'XColor', 'k', 'YColor', 'k', ...
        'Color', 'none');  % Ensure transparent axes background
    axis tight;
    grid on;
end