Avoid root switching using roots

Hello, I am having following problem.

It is a very simple.

First formulate symbolically some matrix .

Then I wish to study . In particular, I am interested only in zeros crossings of imaginary axis while making some sweeps of variables in the matrix, namely

gm2
gm3

So far so simple.

So I get numerator and denumerator of the resulted determinant in each iteration and then find their roots using function roots.

But the problem is, that in the last part of the code

    for pp = 1:size(real_parts_of_poles,1)
        polecross_index = zci(real_parts_of_poles(pp,:));
        
        if ~isempty(polecross_index)
        omega = imag_parts_of_poles(pp,polecross_index);
        MM = [omega, gm2, gm1vec(polecross_index), polecross_index]
        end
    end

when I am detecting these crossing events, I cant proceed like this because obviously roots and their order is mixed during almost each call of roots function.

So instead of smooth evolution like in the subplot below

A have these trajectories from which I cant detect desired crossing events.

Q: is there an alternative how to follow concrete trajectory of pole/zero?

R = 1;
L = 2;
C = 2;

gm1 = 0;
gm2 = 0;
gm3 = 0;
syms s

gm1vec =0:0.1:2.5;
gm2vec = 0:0.5:2;

Z1 = nan(length(gm2vec),length(gm1vec))

real_parts_of_poles = [];
imag_parts_of_poles = [];


for m = 1:length(gm2vec)
    for n = 1:length(gm1vec)
        gm2 = gm2vec(m);
        gm3 = gm1vec(n);
        Y = [1/R + 2/(s*L), -2/(s*L), 0,0,0,0,0,0;
            -2/(s*L), 2/(s*L) + s*C + 1/(s*L), -1/(s*L), 0,0,0,gm3,0;
            0, -1/(s*L), s*C + 2/(s*L), -1/(s*L), 0,gm2,0,0;
            0,0, -1/(s*L), s*C + 2/(s*L), -1/(s*L) + gm1, 0,0,0;
            0,0,0, -1/(s*L), s*C + 2/(s*L), -1/(s*L), 0,0;
            0,0,0,0,-1/(s*L),s*C + 2/(s*L), -1/(s*L), 0;
            0,0,0,0,0, -1/(s*L),s*C + 2/(s*L), -1/(s*L);
            0,0,0,0,0,0, -2/(s*L), 1/R + 2/(s*L);
            ];
        
        dY = det(inv(Y)); 
        
         [N,D] = numden(dY);
        % resultant(N,D,s)
        
        sys = tf(sym2poly(N),sym2poly(D));
        P = pole(sys)
        
        %%% SURFACE PLOTS 13.9.2021
        RR = imag(roots(sym2poly(D)));
        rr = real(roots(sym2poly(D)));
        real_parts_of_poles = round([real_parts_of_poles rr],2);
        imag_parts_of_poles = round([imag_parts_of_poles RR],2);
        
        OO = round(abs(RR(rr>=0)),2);
        
        if isempty(OO)
            Z1(m,n) = nan;
        else %&& (length(OO)<2)
            Z1(m,n) = min(OO(1));
%             counter_first_crossing_first_surface = 1;           
            
            %             if ~counter_first_crossing_second_surface && (length(OO)==2)
            %             Z2(m,n) = OO(2);
            %             counter_first_crossing_second_surface = 1
            %             end
            %
            %             if ~counter_first_crossing_third_surface && (length(OO)==3)
            %             Z3(m,n) = OO(3);
            %             counter_first_crossing_third_surface = 1
            %             end
        end
%         
        
        
%         %
        figure(1)
        subplot(2,1,1)
        plot(real(roots(sym2poly(D))),imag(roots(sym2poly(D))),'x','color',ps{m})
        hold on
        plot(real(roots(sym2poly(N))),imag(roots(sym2poly(N))),'o','color',ps{m})
        
        grid on
        xlabel('$\Re (s)$','interpreter','latex')
        ylabel('$\Im (s)$','interpreter','latex')
        % pause(0.002)
        hAxes = gca;
        hAxes.TickLabelInterpreter = 'latex';
        ax = gca;                        % gets the current axes
        ax.XAxisLocation = 'origin';     % sets them to zero
        ax.YAxisLocation = 'origin';     % sets them to zero
        ax.Box = 'off';                  % switches off the surrounding box
        axis equal
%         
%         
%         
        figure(1)
        subplot(2,1,2)
        plot(gm1vec(n),real(roots(sym2poly(D))),'x','color',ps{m})
        hold on
        plot(gm1vec(n),real(roots(sym2poly(N))),'o','color',ps{m})
        ylabel('$\Re \mathrm{det}\,\mathbf{Z}(s)_{\Omega_0/\Omega_\infty}$','interpreter','latex')
        xlabel('$g_{\mathrm{m}3}$','interpreter','latex')
        grid on
        ylim([-0.2, 0.1])
        
    end
    %%
    %real_parts_of_poles;
    %imag_parts_of_poles;
    
    %figure
    %plot(gm1vec,real_parts_of_poles(6,:))
    
    for pp = 1:size(real_parts_of_poles,1)
        polecross_index = zci(real_parts_of_poles(pp,:));
        
        if ~isempty(polecross_index)
        omega = imag_parts_of_poles(pp,polecross_index);
        MM = [omega, gm2, gm1vec(polecross_index), polecross_index]
        end
    end
    
    
end