Getting strange vertical lines when plotting roots of biquadratic equation/ dispersion relation

Hello,

I'm simply trying to plot the real and imaginary parts of the two positive roots of a biquadratic equation/dispersion relation, p^4 + C3*p^2 + C5 = 0, the plot looks correct apart from these weird vertical lines which seems to be from the roots changing sign...I've racked my brain and I am not sure of what I'm doing wrong. I've put in a divide by zero check, but it doesn't appear that dividing by zeros is the problem.

Here is the plot I'm getting:

and here is what it's supposed to look like:

Thank you in advance for any input!!

%% constants %%
%absolute value of electron charge  
e = 1.60217662E-19; %[C]
%speed of light
c = 299792458; %[m/s]
boltzmann = 1.380649E-23; %[J/K]
m_e = 9.0194E-31; %[kg]
m_i = 6.6335209E-26; %[kg] %argon
% permittivity of free space 
eps_0 = 8.8542E-12; %[F/m]
steps = 1e+3;
B = linspace(0.0002,0.1,steps);%[TESLA]
e_coll_freq = 9e+4;
ion_coll_freq = 6e+5;
dens = 5e+18; %[m^-3]
freq = 8e+6;
w = 2*pi*freq; 
L = 0.8; %[m]
n = 2;
k_parallel = n*pi/L;
N_para = k_parallel*c/w;
% electron cyclotron frequency
w_ce = e.*B/m_e; %[Hz]
% ion cyclotron frequency
w_ci = e.*B/m_i; %[Hz]
% plasma electron frequency
w_pe = sqrt((dens*e^2)/(eps_0*m_e));
% plasma ion frequency 
w_pi = sqrt((dens*e^2)/(eps_0*m_i));
w_LH = sqrt(((w_ce.*w_ci).^2 + (w_pi.*w_ce).^2)./(w_pe^2 + w_ce.^2));
%%%%%

%%% check for dividing by zero %%%
e_denom = (w_ce.^2 - ((1 + 1i*e_coll_freq/w)^2)*w^2);
divZeroElements = find(abs(e_denom)<0.01);
e_denom(divZeroElements) = NaN;
i_denom = (w_ci.^2 - ((1 + 1i*ion_coll_freq/w)^2)*w^2);
divZeroElements2 = find(abs(i_denom)<0.01); 

eps1 = (1 + ((1 + 1i*e_coll_freq/w)*w_pe^2)./e_denom + ((1 + 1i*ion_coll_freq/w)*w_pi^2)./i_denom);

eps2 = (-((w_pe^2).*(w_ce./w))./e_denom + ((w_pi^2)*(w_ci./w))./i_denom);

eps3 = (1 - w_pe^2/(((1 + 1i*e_coll_freq/w)^2)*w^2) - w_pi^2/(((1 + 1i*ion_coll_freq/w)^2)*w^2));

Alpha = eps1 - N_para^2 - (eps2.^2)./eps1;

Beta = eps3.*(1-((N_para^2)./eps1));

Delta = N_para.*eps2./eps1;

Gamma = (N_para.*eps2.*eps3)./eps1;

%coefficients of 4th order polynomial 
C1 = 1;
C2 = 0;
C3 = (-1).*(Alpha+Beta);
C4 = 0;
C5 = (Alpha.*Beta) - (Gamma.*Delta);

ROOTS1(:,1) = sqrt((-C3 + sqrt(C3.^2 - 4.*C5))./2);
ROOTS2(:,1) = sqrt((-C3 - sqrt(C3.^2 - 4.*C5))./2);
ROOTS3(:,1) = -sqrt((-C3 + sqrt(C3.^2 - 4.*C5))./2);
ROOTS4(:,1) = -sqrt((-C3 - sqrt(C3.^2 - 4.*C5))./2);

% normalized x_axis
x_axis = w_LH./w;

figure()
plot(x_axis,imag(ROOTS1),'--k')
hold on
plot(x_axis,real(ROOTS1))
plot(x_axis,real(ROOTS2))
plot(x_axis,imag(ROOTS2),'--b')
ylim([-1.5e+4 1.5e+4])
xlabel('\omega_{LH}/\omega')
ylabel('p')

Answer 1

[Katey Stevenson]

I fixed it, thankfully, using rootshuffle.m https://www.mathworks.com/matlabcentral/fileexchange/29464-rootshuffle-m?s_tid=srchtitle !!! Apparently, it is because matlab finds the roots, but doesn't necessarily put them in order, thereby causing sign flipping of the roots!!!