Hello Shiv,
I would proceed as follows. The figure is not similar, but it may be useful to see the procedure.
Best,
Alberto
% Definitions
% * Constants
C_t = 1;
C_cr = 1;
v0 = 0;
va = 1;
w_j = 0;
% * Domain
k = linspace(0, 0.7);
% Compute polynomial coefficients
a0 = va^4 * k.^4 .* (k.^2 * (C_t^2 + C_cr^2) - w_j^2);
a1 = 0;
a2 = (va^2 * k.^4 + (k.^2 * (C_t^2 + C_cr^2) - w_j^2) .* (2 * va^2 * k.^2 + 4 * v0^2 * k.^4));
a3 = 0;
a4 = k.^2 .* (4 * (va^2 + 2 * v0^2 .* k.^2) + C_t^2 + C_cr^2);
a5 = 0;
a6 = 1;
% Obtain w
for i = length(k):-1:2
p = [a0(i), a1, a2(i), a3, a4(i), a5, a6];
w(:, i) = roots(p);
end
% Plot the real part of all solutions of w
plot(k, real(w));
