MATLAB: Vpasolve for three nonlinear equations inside a for loop

for loopnon-linear equationsvpasolve

I am trying to solve the system of the three non-linear equations by using vpasolve for different values of parameter U that varies from 0-0.5, while other parameters are fixed. My code is just showing output for U=0, then it's showing an error(Second argument must be a vector of symbolic variables).
clc
clear all;
syms x y w z eb
t = 0.2./pi;
d = 0.2./pi;  
U = 0;
 while (U<0.5)
e=U./2;
a = U./(pi.*t);
   f = imag((U./d).*((((t./d)./(sqrt(1-z.^2)))+w)./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
       sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2)));
   g = imag((z+(( -e + U.*x)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
       sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2));
   h = imag((z+((-e + U.*y)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
       sqrt(1-z.^2))-(t./d).^2-((-e + U.*y)./d).^2-w.^2)); 
   
        s=int(f,z,-Inf,0);
        u=int(g,z,-Inf,0);
        v=int(h,z,-Inf,0); 
        
       eq1=w-(-1./pi).*s==0;
       eq2=y-(-1./pi).*u==0;
       eq3=x-(-1./pi).*v==0;
         
% sol = vpasolve(eqs,vars);   
[x,y,w] = vpasolve([eq1, eq2, eq3],[x,y,w],[1 0 0]);
%[sol.x sol.y sol.w]
%solutions = [solx,soly,solw]
n_up = double(x);
n_down = double(y);
d_ind = double(w);
m = abs(n_up-n_down);
eqn=((eb.*(1+2.*(t./d)./sqrt(1-eb)))-(2.*(t./d).*d_ind./sqrt(1-eb))-(t./d).^2-(( -e + U.*n_up)./d).^2-d_ind.^2);
e_abs = vpasolve(eqn,eb);
e_abs1=sqrt(e_abs);
e_abs2=-sqrt(e_abs);
weight = ((1./2).*(1-(e_abs1).^2).*(((sqrt(1-e_abs1.^2)).*(1+((( -e + U.*n_up)./d)./(e_abs1))))./(((1-(e_abs1).^2)...
    .*((sqrt(1-e_abs1.^2))+2.*(t./d)))+((t./d).*e_abs1.^2)+(t./d).*(d_ind./d))));
fprintf('%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f  %8.4f %8.4f\n', [a,n_up,n_down,m,d_ind,e_abs1,e_abs2,weight]');
if (U==0.5)
     break
end
U = U+0.1;
 end

Best Answer

In this line
[x,y,w] = vpasolve([eq1, eq2, eq3],[x,y,w],[1 0 0]);
you are overwriting the values of x, y, and w and converting them from symbolic to numeric. Use different variable names. See the following code
clc
clear all;
syms x y w z eb
t = 0.2./pi;
d = 0.2./pi;
U = 0;
while (U<0.5)
    e=U./2;
    a = U./(pi.*t);
    f = imag((U./d).*((((t./d)./(sqrt(1-z.^2)))+w)./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
        sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2)));
    g = imag((z+(( -e + U.*x)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
        sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2));
    h = imag((z+((-e + U.*y)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
        sqrt(1-z.^2))-(t./d).^2-((-e + U.*y)./d).^2-w.^2));
    
    s=int(f,z,-Inf,0);
    u=int(g,z,-Inf,0);
    v=int(h,z,-Inf,0);
    
    eq1=w-(-1./pi).*s==0;
    eq2=y-(-1./pi).*u==0;
    eq3=x-(-1./pi).*v==0;
    
    % sol = vpasolve(eqs,vars);
    [xv,yv,wv] = vpasolve([eq1, eq2, eq3],[x,y,w],[1 0 0]);
    %[sol.x sol.y sol.w]
    %solutions = [solx,soly,solw]
    n_up = double(xv);
    n_down = double(yv);
    d_ind = double(wv);
    m = abs(n_up-n_down);
    eqn=((eb.*(1+2.*(t./d)./sqrt(1-eb)))-(2.*(t./d).*d_ind./sqrt(1-eb))-(t./d).^2-(( -e + U.*n_up)./d).^2-d_ind.^2);
    e_abs = vpasolve(eqn,eb);
    e_abs1=sqrt(e_abs);
    e_abs2=-sqrt(e_abs);
    weight = ((1./2).*(1-(e_abs1).^2).*(((sqrt(1-e_abs1.^2)).*(1+((( -e + U.*n_up)./d)./(e_abs1))))./(((1-(e_abs1).^2)...
        .*((sqrt(1-e_abs1.^2))+2.*(t./d)))+((t./d).*e_abs1.^2)+(t./d).*(d_ind./d))));
    fprintf('%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f  %8.4f %8.4f\n', [a,n_up,n_down,m,d_ind,e_abs1,e_abs2,weight]');
    if (U==0.5)
        break
    end
    U = U+0.1;
end
Related Question