Solving multiple cubic equations in loop

I have 'l' cubic equations, and I want to extract the imaginary part of the root of these equations.

Below is my code. The part in the for loop is where i'm having trouble. I don't understand how to collect the imaginary part of the roots.

close all
clear all
clc

global a0 a1 a2 a3 We U


rhog = 1.2;
rhol = 1.0e3;
r = rhog/rhol;
sigma = 0.072; % surface tension
D = 0.0025; % diameter of jet
U = 1500;    %cm/s
q = 3.886*10e-4 ;%cm2
x=1/2;    %cm
h = q/x; %cm
% non-dimensional
We = 1/(rhol*(U/100)^2*h/(100*sigma)) % 1/Weber number

%------End of User defined constants------------------4

omega_max = 20000;
for l= 1:omega_max
   w(l) = l+1;        %frequency
   a0(l) = (-r*w(l)^2/((U^2)*h*(1-We)));   %Equation: a3*x^3 + a2*x^2 + a1*x + a0 = 0
   a1(l) = (-w(l)^2/((U^2)*(1-We)));
   a2(l) = (2*w(l)/(U*(1-We)));
   a3(l) = 1;
   k(l) = roots([a3(l) a2(l) a1(l) a0(l)]);
   k_im(l) = imag(k(l))*ij;
end 

plot( w/(2*22/7),k_im,'r*')  
xlabel('frequency(Hz)')
ylabel('k imaginary (cm^-^1)')

[Vladimir Sovkov]

A cubic equation has 3 complex roots provided by the function "roots" but you are trying to deliver it to one scalar element of a vector k(i), that is erroneous.

Answer 1

[Walter Roberson]

   k(l) = roots([a3(l) a2(l) a1(l) a0(l)]);

That will fail because there are going to be 3 roots of a cubic and you are trying to write all 3 into a single slot k(l)

k(l, :) = roots([a3(l) a2(l) a1(l) a0(l)]);
k_im(l, :) = imag(k(l,:))*ij;

[Soumya Kedia]

Yes thanks!