One way is to collect all the points to be plotted into two matrices, one for x-coordinates and one for y-coordinates, remove the imaginary and non-positive y values (not just set them to NaN, since that creates gaps), sort the resulting vectors according to x, and plot a single line.
clear all ;
syms V omega
xi = 0.2;
gamma = 0.3;
chi = 0.1e-1;
lambda =1;
eta = 1;
alpha =0.7;
beta = 0.1;
% omega=1.15;
F=0.3;
% Given parameters
k1 = -(1/2).*xi;
k2 = chi./(2.*omega);
k3 = F./(2.*omega);
k4 = (-omega.^2+1)./(2.*omega);
k5 = 3.*gamma./(8.*omega);
k6 = (1./2).*alpha-(1./2).*lambda;
k7 = -(1./8).*beta;
k8 = (1./2).*omega;
k9 = (eta.^2-omega.^2)./(2.*omega);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A1 = -k7.^2;
A2 = -2.*k6.*k7;
A3 = -k6.^2-k9.^2;
A4 = k8.^2;
A5 = k5.^2.*k8.^2;
A6 = 2.*k4.*k5.*k8.^2;
A7 = -2.*k2.*k5.*k8.*k9;
A8 = k1.^2.*k8.^2+k4.^2.*k8.^2;
A9 = 2.*k1.*k2.*k7.*k8;
A10 = 2.*k1.*k2.*k6.*k8-2.*k2.*k4.*k8.*k9;
A11 = -k3.^2.*k8.^2;
A12 = k2.^2.*k7.^2;
A13 = 2.*k2.^2.*k6.*k7;
A14 = k2.^2.*k6.^2+k2.^2.*k9.^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the system of equations
EQ = V^22*A1^4*A5+4*V^20*A1^3*A2*A5+(4*A1^3*A3*A5+6*A1^2*A2^2*A5)*V^18+(-A1^3*A4*A6+12*A1^2*A2*A3*A5+4*A1*A2^3*A5)*V^16+(-3*A1^2*A2*A4*A6+6*A1^2*A3^2*A5+12*A1*A2^2*A3*A5+A2^4*A5)*V^14+(-3*A1^2*A3*A4*A6+A1^2*A4^2*A7-3*A1*A2^2*A4*A6+12*A1*A2*A3^2*A5+4*A2^3*A3*A5)*V^12+(A1^2*A4^2*A8-6*A1*A2*A3*A4*A6+2*A1*A2*A4^2*A7+4*A1*A3^3*A5-A2^3*A4*A6+6*A2^2*A3^2*A5)*V^10+(2*A1*A2*A4^2*A8-3*A1*A3^2*A4*A6+2*A1*A3*A4^2*A7-A1*A4^3*A9-3*A2^2*A3*A4*A6+A2^2*A4^2*A7+4*A2*A3^3*A5)*V^8+(2*A1*A3*A4^2*A8-A1*A4^3*A10+A2^2*A4^2*A8-3*A2*A3^2*A4*A6+2*A2*A3*A4^2*A7-A2*A4^3*A9+A3^4*A5+A4^4*A12)*V^6+(-A1*A4^3*A11+2*A2*A3*A4^2*A8-A2*A4^3*A10-A3^3*A4*A6+A3^2*A4^2*A7-A3*A4^3*A9+A4^4*A13)*V^4+(-A2*A4^3*A11+A3^2*A4^2*A8-A3*A4^3*A10+A4^4*A14)*V^2-A3*A4^3*A11;
% Solve the system of equations
sol = solve(EQ,V);
% % Display the solutions
% disp('Solutions for V:');
% disp(sol);
% Plot the solutions versus omega
omega_values = linspace(0.1, 1.8, 300); % adjust the range accordingly
figure(11);
hold on;
% collect x and y points for all solutions in matrices,
% one column per solution
Nsol = numel(sol);
x_plot = repmat(omega_values(:),1,Nsol);
y_plot = zeros(numel(omega_values),Nsol);
for i = 1:Nsol
y_plot(:,i) = double(subs(sol(i), omega, omega_values));
end
% remove imaginary and non-positive values
% (makes x_plot and y_plot into row vectors)
idx = imag(y_plot) ~= 0 | y_plot <= 0;
x_plot(idx) = [];
y_plot(idx) = [];
% sort by increasing x_plot values
[x_plot,idx] = sort(x_plot);
y_plot = y_plot(idx);
% plot
plot(x_plot, y_plot, 'b-', "linewidth", 2);
hold off;
box on;
set(gca, "linewidth", 1.5, 'FontSize', 14, 'FontWeight', "bold");
xlabel('\omega', 'fontname', 'Times New Roman', 'FontSize', 16);
ylabel('V', 'fontname', 'Times New Roman', 'FontSize', 16);
