function plotPhaseDiagramWithBoundaries(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, PhaseBoundary, TitleString)                                         

    % plotPhaseDiagramWithBoundaries - Plots a phase diagram with overlaid interpolated boundaries
    %
    % Inputs:
    %   PhaseDiagramMatrix     - 2D matrix of phase values
    %   NumberOfAtoms          - Vector of atom numbers (x-axis)
    %   ScatteringLengths      - Vector of scattering lengths (y-axis)
    %   PhaseBoundary_Untilted - Object containing selected boundary points
    %   TitleString            - Title string for the plot
    %
    % This function displays the interpolated phase diagram and overlays
    % boundary curves extracted from user-selected points.

    % Extract boundary point sets
    rawPoints = PhaseBoundary.SelectedPoints.values;
    if iscell(rawPoints)
        pointSets = rawPoints;
    else
        pointSets = {rawPoints};
    end

    regionNames = ["Unmodulated", "SSD", "Stripe", "Labyrinth", "Honeycomb"];
    regionColors = [
        0.8 0.8 0.8;
        0.2 0.6 1.0;
        0.2 0.8 0.2;
        1.0 0.6 0.2;
        0.8 0.2 0.8
    ];

    [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');

    fig = figure('Color', 'w'); clf;
    set(fig, 'Position', [100, 100, 1050, 800], 'Renderer', 'opengl');
    ax = axes('Parent', fig, 'Position', [0.1 0.15 0.75 0.8]);

    imagesc(ax, [min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST)], ...
                [min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE)], Mq);
    set(ax, 'YDir', 'normal');
    colormap(ax, regionColors);
    cb = colorbar('Ticks', 0:4, 'TickLabels', regionNames, 'FontSize', 12);
    cb.Color = 'k';  % set colorbar text color
    clim(ax, [0 4]);
    hold(ax, 'on');
    xlabel(ax, 'Number of Atoms', 'FontSize', 16);
    ylabel(ax, 'Scattering Length (\times a_o)', 'FontSize', 16,'Interpreter', 'tex');
    t = title(ax, TitleString, 'FontSize', 18, 'Interpreter', 'tex');
    t.Color = 'k';  % force title color to black
    set(ax, 'FontSize', 16, 'Color', 'none');
    axis(ax, 'tight');
    grid(ax, 'on');

    % Overlay curves
    hold on;
    for i = 1:numel(pointSets)
        pts = pointSets{i};
        x = pts(:,1) * 1e5;  % convert back to # of atoms scale
        y = pts(:,2);

        % Interpolation
        [x, idx] = sort(x); y = y(idx);
        dist = sqrt(diff(x).^2 + diff(y).^2);
        t = [0; cumsum(dist)];
        t = t / t(end);
        t_fine = linspace(0, 1, 1000);
        x_smooth = interp1(t, x, t_fine, 'makima');
        y_smooth = interp1(t, y, t_fine, 'makima');

        % Plot smooth curve
        plot(x_smooth, y_smooth, 'k-', 'LineWidth', 2);
    end

    set(gca, 'FontSize', 16, ...
            'XColor', 'k', 'YColor', 'k', ...
            'Color', 'none');
    axis tight;
    grid on;

end