Error with vpasolve - arithmetical expression expected

Hello!

I have to solve a symbolic equation in 2 variables (the real and imaginary parts of a complex variable) inside of a for loop.

I have already tested the script with another equation, so I am now adapting it to solve this new equation from a new model, hence only changing the parameters and equation and leaving everything else unchanged.

However, now the solver is not working and I receive the following error message:

clear
close all
clc

%Constant parameters%

a = 0.005; %sheet half thickness%
eps0 = 8.9*10^-12; %vacuum dielectric constant%
eps_2 = 80; %water relative dielectric constant% 
epsg = 1.006; %air relative dielectric constant%
epsilon_G = eps0*epsg; %air dielectric constant%
epsilon_L = eps0*eps_2; %water dielectric constant%
gamma = 0.0728; %water-air surface tension%
mu = 8.94*10^-4; %water viscosity%
rho_2 = 1000; %water density%
rhog = 1.2250; %air density%
sigma = 500; %water conductivity%

%Other parameters and variables%
U = 1; %liquid velocity%

%Plotted conditions%

We = 400; %Weber number%
Re = 1000; %Reynolds number%
rho = 0.0012; %density ratio%
tau = 0.01; %electrical relaxation time/hydrodynamic time%
Eu  = 2; %Euler number%
eps = 0.0125; %epsilon_G/epsilon_L%
D = 40; %dimensionless distance d/a%


%Dispersion relation part%

syms omega_r omega_i

k_ = 1.1; %ending point for k%
delta = 0.01; 
sizek = numel(0:delta:k_); %dimension of the k vector%

omega_r_arr = zeros(sizek,1);
omega_i_arr = zeros(sizek,1);
k = zeros(sizek,1);
k(1) = 0;


for j = 1:sizek
    
    
    omega = omega_r+1i*omega_i; %complex number omega%
    l_2 = k(j)^2+(rho_2/mu)*(omega+U*1i*k(j)); %l^2%
    l = sqrt(l_2); %l%
    
    V0_d_D32 = -((l_2-k(j)^2)/(k(j)^2-l_2))*2*coth(k(j))*cosh(k(j))*(k(j)*Eu)/D^2 ...
        +((4*k(j)^2*l)/(k(j)^2-l_2))*(cosh(k(j))/tanh(l))*(Eu/D^2) ...
        -2*sinh(k(j))*(k(j)*Eu/D^2)
    
    D33 = -((l_2+k(j)^2)^2/k(j))*tanh(k(j))/Re^2 ...
        +(4*k(j)^2*l/Re^2)*tanh(l)...
        -(k(j)^2/We) ...
        -(rho*omega^2/k(j)) ...
        -((l_2+k(j)^2)^2/(k(j)^2-l_2))*(k(j)*Eu*coth(k(j))/D^2) ...
        +(2*k(j)^2*l/(k(j)^2-l_2))*(Eu*coth(l)/D^2) ...
        +Eu/D^3;
    
    
    DR = (tau-(l_2-k(j)^2)/Re)*(-2*cosh(k(j)))*(-((l_2+k(j)^2)^2/k(j))*tanh(k(j))/Re^2 ...
        +(4*k(j)^2*l/Re^2)*tanh(l)...
        -(k(j)^2/We) ...
        -(rho*omega^2/k(j)) ...
        -((l_2+k(j)^2)^2/(k(j)^2-l_2))*(k(j)*Eu*coth(k(j))/D^2) ...
        +(2*k(j)^2*l/(k(j)^2-l_2))*(Eu*coth(l)/D^2) ...
        +Eu/D^3) ...
        +eps*((l_2-k(j)^2)/Re)*(cosh(k(j))/(k(j)*D*sinh(k(j))))*(-((l_2-k(j)^2)/(k(j)^2-l_2))*2*coth(k(j))*cosh(k(j))*(k(j)*Eu)/D^2 ...
        +((4*k(j)^2*l)/(k(j)^2-l_2))*(cosh(k(j))/tanh(l))*(Eu/D^2) ...
        -2*sinh(k(j))*(k(j)*Eu/D^2)) ...
        +(tau-(l_2-k(j)^2)/Re)*(-((l_2-k(j)^2)/(k(j)^2-l_2))*2*coth(k(j))*cosh(k(j))*(k(j)*Eu)/D^2 ...
        +((4*k(j)^2*l)/(k(j)^2-l_2))*(cosh(k(j))/tanh(l))*(Eu/D^2) ...
        -2*sinh(k(j))*(k(j)*Eu/D^2)) ...
        -eps*((l_2-k(j)^2)/Re)*(2*cosh(k(j)))*(-((l_2+k(j)^2)^2/k(j))*tanh(k(j))/Re^2 ...
        +(4*k(j)^2*l/Re^2)*tanh(l)...
        -(k(j)^2/We) ...
        -(rho*omega^2/k(j)) ...
        -((l_2+k(j)^2)^2/(k(j)^2-l_2))*(k(j)*Eu*coth(k(j))/D^2) ...
        +(2*k(j)^2*l/(k(j)^2-l_2))*(Eu*coth(l)/D^2) ...
        +Eu/D^3);
    
    
    [wr(j), wi(j)] = vpasolve([real(DR)==0,imag(DR)==0],[omega_r, omega_i], [3*10^-3 3*10^-3])
    
    
    omega_r_arr(j) = wr(j);
    omega_i_arr(j) = wi(j);
    
    k(j+1) = k(j)+delta;
    
end

V0_d_D32 = 

NaN

Error using mupadengine/feval_internal
Arithmetical expression expected.

Error in sym/vpasolve (line 172)
sol = eng.feval_internal('symobj::vpasolve',eqns,vars,X0);

k=0:delta:k_;
plot(k,omega_i_arr)
grid on

I don't understand what it means and how to correct it, can anybody please help me? I am attaching the matlab code to the question.

[Torsten]

See above.

V0_d_D32 becomes NaN.

[RITA PIA SESSA]

Ok thank you! I might have copied something wrong from the source then, thank you again for your answer