Refinement of your code:
syms a omega
% Given parameters
omega_n =1;
mu = 0.1;
alpha = 0.2;
chi = 0.4;
beta = 0.1;
kappa = 0.15;
F1 = 0.3;
F3 =0;
sigma = omega - omega_n;
B = -(1/8)*F3/omega_n^2;
H1 = -4*beta^2*F1 - 4*F1*omega_n^2;
H2 = 6*B*alpha*beta^2 + 6*B*alpha*omega_n^2;
H3 = 4*beta^2*mu*omega_n + 4*beta*chi*kappa*omega_n + 4*mu*omega_n^3;
H4 = 3*alpha*beta^2 + 3*alpha*omega_n^2;
H5 = -24*B^2*alpha*beta^2 - 24*B^2*alpha*omega_n^2 + 8*beta^2*sigma*omega_n - 4*chi*kappa*omega_n^2 + 8*sigma*omega_n^3;
G1 = H4^6;
G2 = -3*H2^2*H4^4 - 6*H4^5*H5;
G3 = 3*H2^4*H4^2 + 12*H2^2*H4^3*H5 + 3*H3^2*H4^4 + 15*H4^4*H5^2;
G4 = -H1^2*H4^4 - H2^6 - 6*H2^4*H4*H5 - 6*H2^2*H3^2*H4^2 - 18*H2^2*H4^2*H5^2 - 12*H3^2*H4^3*H5 - 20*H4^3*H5^3;
G5 = -H1^2*H2^2*H4^2 + 4*H1^2*H4^3*H5 + 3*H2^4*H3^2 + 3*H2^4*H5^2 + 12*H2^2*H3^2*H4*H5 + 12*H2^2*H4*H5^3 + 3*H3^4*H4^2 + 18*H3^2*H4^2*H5^2 + 15*H4^2*H5^4;
G6 = -2*H1^3*H2*H4^2 + 2*H1^2*H2^4 + 2*H1^2*H2^2*H4*H5 - 2*H1^2*H3^2*H4^2 - 6*H1^2*H4^2*H5^2 - 3*H2^2*H3^4 - 6*H2^2*H3^2*H5^2 - 3*H2^2*H5^4 - 6*H3^4*H4*H5 - 12*H3^2*H4*H5^3 - 6*H4*H5^5;
G7 = 4*H1^3*H2*H4*H5 - H1^2*H2^2*H3^2 - H1^2*H2^2*H5^2 + 4*H1^2*H3^2*H4*H5 + 4*H1^2*H4*H5^3 + H3^6 + 3*H3^4*H5^2 + 3*H3^2*H5^4 + H5^6;
G8 = -H1^4*H2^2 + 2*H1^3*H2*H3^2 - 2*H1^3*H2*H5^2 - H1^2*H3^4 - 2*H1^2*H3^2*H5^2 - H1^2*H5^4;
% Define the system of equations
EQ = a^14*G1 + a^12*G2 + a^10*G3 + a^8*G4 + a^6*G5 + a^4*G6 + a^2*G7 + G8;
% Solve the system of equations
sol = solve(EQ, a);
% Display the solutions
disp('Solutions for a:');
disp(sol);
% Plot the solutions versus omega
omega_values = linspace(0,1.8,300); % adjust the range accordingly
figure(2);
hold on;
for i = 1:length(sol)
a_values = double(subs(sol(i), omega, omega_values));
a_values(imag(a_values)~=0) = nan;
plot(omega_values, a_values, 'DisplayName', ['Solution ' num2str(i)]);
end
hold off;
xlabel('\omega');
ylabel('a');
title('Solutions versus \omega');
legend('show');
