how to solve this error ?

In this code, MATLAB gives me the following issue: For the solving function, how can I correct this error?

Error using sym/vpasolve:

Unable to find variables in equations.

Error in ni_water_simulationburgemanforchatgbt (line 99):

`solutions = vpasolve(eq, epsilon_mg, min(epsi, epsh), max(epsi, epsh));`

clear all
lambda = 500*10^-3:10*10^-3:2500*10^-3; % in microns
A1= 1.4182;
B1= 0.021304;
n_pc=sqrt(1+(A1.*lambda.^2)./(lambda.^2-B1));
wp= 15.92; % eV
f0 = 0.096;
Gamma0 = 0.048; % eV
f1 = 0.100;
Gamma1 = 4.511; % eV
omega1 = 0.174; % eV
f2 = 0.135;
Gamma2 = 1.334; % eV
omega2 = 0.582; % eV
f3 = 0.106;
Gamma3 = 2.178; % eV
omega3 = 1.597; % eV
f4 = 0.729;
Gamma4 = 6.292; % eV
omega4 = 6.089; % eV
OmegaP = sqrt(f0) * wp; % eV
eV = 4.13566733e-1 * 2.99792458 ./ lambda;
epsilon = 1 - OmegaP^2 ./ (eV .* (eV + 1i * Gamma0));
epsilon = epsilon + f1 * wp^2 ./ ((omega1^2 - eV.^2) - 1i * eV * Gamma1);
epsilon = epsilon + f2 * wp^2 ./ ((omega2^2 - eV.^2) - 1i * eV * Gamma2);
epsilon = epsilon + f3 * wp^2 ./ ((omega3^2 - eV.^2) - 1i * eV * Gamma3);
epsilon = epsilon + f4 * wp^2 ./ ((omega4^2 - eV.^2) - 1i * eV * Gamma4);
n = real(sqrt(epsilon));
k = imag(sqrt(epsilon));
nMetal = n+1i.*k;
%NWATER 
NWATER = 0.0738.*lambda.^6 - 0.6168.*lambda.^5 + 2.0263.*lambda.^4 - 3.3315.*lambda.^3 + 2.8708.*lambda.^2 - 1.2367.*lambda + 1.5411;
nAir=1+0.05792105./(238.0185-(lambda.^(-2)))+0.00167917./(57.362-(lambda.^(-2)));
% permetivity 
epsMetal = nMetal.^2;
epsWATER=NWATER.^2; 
%parmenter for each axes
%1 specific volume 
ni=n_pc;
nt=NWATER;
%angels 
thetaI =0;               % incidence angle
thetaI = thetaI/180*pi;
beta=20;
beta = beta/180*pi; 
fi=0;
fi=fi/180*pi;
%%%%%
a=20*10.^-9;
b=5*10.^-9;
d=2*a*cos(beta);
fi=0.3;
fh=1-fi;
Rpp = [];
Rsp = [];
Rps = [];
Rss = [];
Tpp = [];
Tsp = [];
Tps = [];
Tss = [];
% epsilon for compost materals SHAPE Factor 
e= sqrt(1-(b/a)^2);
la=((1-e^.2)/e^2)*((log((1+e)/(1-e)))/2*e^2  -1);
lb=0.5*(1-la);
% Symbolic variable for epsilon_mg
syms epsilon_mg
%%%%%%%%
for i=1:size(lambda,2)
    %snel law §
    thetaT=asin(ni(i)*sin(thetaI)/nt(i));
    k0=2*pi/(lambda(i)*10^-6);
    VX=ni(i)*sin(thetaI);
    epsi=epsMetal(i);
    epsh=epsWATER(i);
    
    % Bruggeman mixing formula
    eq = fi * (epsi - epsilon_mg) / (epsilon_mg + la * (epsi - epsilon_mg)) +(1 - fi) * (epsh - epsilon_mg) / (epsilon_mg + lb * (epsh - epsilon_mg)) == 0;
    % Solve equation with positive imaginary part
    solutions = vpasolve(eq, epsilon_mg, [min(epsi, epsh), max(epsi, epsh)]);
    positive_imaginary = solutions(imag(solutions) > 0);
    if ~isempty(positive_imaginary)
        epsilon_mg = positive_imaginary(1);
    else
        error('No solution with positive imaginary part found at index %d', i);
        
    end
    eps1=epsilon_mg;
    eps2=epsilon_mg;
    eps3=epsilon_mg;
    exx=eps2 +(eps1.*cos(beta).*cos(beta) +eps3.*sin(beta).*sin(beta) -eps2).*cos(fi).*cos(fi);
    exy=0.5.*(eps1.*cos(beta).*cos(beta)  +eps3.*sin(beta).*sin(beta) - eps2).*sin(2.*fi);
    exz=0.5.*(eps3-eps1).*sin(2.*beta).*cos(fi);
    eyy=eps2 +(eps1.*cos(beta).*cos(beta)  +eps3.*sin(beta).*sin(beta) - eps2).*sin(fi).*sin(fi);
    eyz=0.5.*(eps3-eps1).*sin(2.*beta).*sin(fi);
    ezz=eps1 +(eps3 -eps1).*cos(beta).*cos(beta);
    eyx=exy;
    ezx=exz;
    ezy=eyz;
    del = [-VX.*ezx./ezz  1-(VX.*VX)./ezz  -VX.*ezy./ezz 0; exx-(exz.*ezx)./ezz  -VX.*exz./ezz  exy-exz.*ezy./ezz 0; 0 0 0 1; eyx-(eyz.*ezx)./ezz  -VX.*eyz./ezz eyy-(VX.*VX)-(eyz.*ezy)./ezz  0];
    P=exp(1)^(1i.*k0.*d.*del);
    a1 = ni(i).*(nt(i).*P(1,2)-cos(thetaT).*P(2,2))+cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
    a2 = ni(i).*(nt(i).*P(1,2)-cos(thetaT)*P(2,2))-cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
    a3 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))+ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
    a4 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))-ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
    a5 = ni(i).*(nt(i).*cos(thetaT).*P(3,2 )-P(4,2)) +cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
    a6 = ni(i).*(nt(i).*cos(thetaT).*P(3,2) -P(4,2))-cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
    a7 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))+ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
    a8 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))-ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
    b1 = (ni(i).*P(2,2)+cos(thetaI).*P(2,1))./nt(i);
    b2 = (ni(i).*P(2,2)-cos(thetaI).*P(2,1))./nt(i);
    b3 = (P(2,3)-ni(i).*cos(thetaI).*P(2,4))./nt(i);
    b4 = (P(2,3)+ni(i).*cos(thetaI).*P(2,4))./nt(i);
    b5 = ni(i).*P(3,2)+cos(thetaI).*P(3,1);
    b6 = ni(i).*P(3,2)-cos(thetaI).*P(3,1);
    b7 = P(3,3)-ni(i).*cos(thetaI).*P(3,4);
    b8 = P(3,3)+ni(i).*cos(thetaI).*P(3,4);
    rpp=(a1.*a8-a4.*a5)./(a4.*a6-a2.*a8);
    rps=(a3.*a8-a4.*a7)./(a4.*a6-a2.*a8);
    rsp=(a2.*a5-a1.*a6)./(a4.*a6-a2.*a8);
    rss=(a2.*a7-a6.*a3)./(a4.*a6-a2.*a8);
    tpp=b1+b2.*rpp+b3.*rsp;
    tps=b4+b2.*rps+b3.*rss;
    tsp=b5+b6.*rpp+b7.*rsp;
    tss=b8+b6.*rps+b7.*rss;
    Rpp = [Rpp (abs(rpp))^2];
    Rsp = [Rsp (abs(rsp))^2];
    Rps = [Rps (abs(rps))^2];
    Rss = [Rss (abs(rss))^2];
    Tpp = [Tpp  (nt(i) .* cos(thetaI) ./ ni(i) .* cos(real(thetaT))).*(abs(tpp )^2)];
    Tss = [Tss  (nt(i) .* cos(real(thetaT)) ./ ni(i) .* cos(thetaI)).*(abs(tss )^2)];
    Tsp = [Tsp  (nt(i) .* cos(real(thetaT)) ./ ni(i) .* cos(thetaI)).*(abs(tsp )^2)];
    Tps = [Tps  (nt(i) .* cos(real(thetaT)) ./ ni(i) .* cos(thetaI)).*(abs(tps )^2)];
    T=(Tpp+Tsp+Tps+Tss)/2;
    R=(Rpp+Rsp+Rps+Rss)/2;
    A=1-R-T;
end 
figure(1)
hold on
plot(lambda*10^3,A*100,'LineWidth', 3)
hold on
xlabel('Wavelength (nm)','fontsize',16)
hold on
ylabel('Absorbtion','fontsize', 16)
figure(2)
hold on
plot(lambda*10^3,T*100,'LineWidth', 3)
hold on
xlabel('Wavelength (nm)','fontsize',16)
hold on
ylabel('Transmission','fontsize', 16)

Answer 1

[Walter Roberson]

    solutions = vpasolve(eq, epsilon_mg, [min(epsi, epsh), max(epsi, epsh)]);
    positive_imaginary = solutions(imag(solutions) > 0);
    if ~isempty(positive_imaginary)
        epsilon_mg = positive_imaginary(1);
    else
        error('No solution with positive imaginary part found at index %d', i);
        
    end

After the first iteration of this, either you will get error() triggered or else epsilon_mg will be assigned the value of the first solution. But the solution is the output of vpasolve() and so is a symbolic expression with no variables in it, because vpasolve() can only return symbolic numbers or empty, never expressions involving symbolic variables.

So, the second iteration of the loop, epsilon_mg is no longer an unresolved symbolic variable and is instead some symbolic number. eq involves calculation using that symbolic number, so eq will be some symbolic number that does not include any symbolic variables. The chances are low that eq and epsilon_mg both happen to be equivalent to 0, so the vpasolve() will fail.

In other words, your

    syms epsilon_mg

should be inside of the loop if you are going to overwrite epsilon_mg the way you are.

[aideen]

it did work ,Thanks for your help in solving this problem.