vpasolve of a equation with four solutions.

Hello, I am working with the fresnell equation (to generate the surface of biaxic materials) whose solutions are four and symetrical. The problem is I only obtain 1 solution in 2 declared graphics and both solutions are equal.

How is it possible to obtain the other solutions?

clear all
syms n

Nx=1.630;
Ny=1.634;
Nz=1.638;

fprintf('Vamos a resolver\n');
%x=sin(phi)*cos(theta)
%y=sin(phi)*sin(theta)
%z=cos(phi)
S=[];
Coord=[];
theta=linspace(0,pi,20);
phi=linspace(-pi/2,pi/2,20);
[X,Y]=meshgrid(theta,phi);
for i=1:length(X)
    fprintf('%d\n',i)
    for j=1:length(Y)
            %S=[sin(Y(j))*cos(X(i)) sin(Y(j))*sin(X(i)) cos(Y(j))];       
            %x=S(1)/norm(S);y=S(2)/norm(S);z=S(3)/norm(S);
            eqn=((sin(Y(j)).*cos(X(i))).^2/(n.^2-Nx^2)+(sin(Y(j)).*sin(X(i))).^2/(n.^2-Ny^2)+(cos(Y(j))).^2./(n.^2-Nz^2)).*(n.^2)==1;
            sol=double(vpasolve(eqn,n,[-Inf,Inf]));
            A=sol>0 & sol<10;
            B=real(sol(A));
                if length(B)==1
                    Z1(i,j)=B;
                    Z2(i,j)=B;
                else
                    Z1(i,j)=B(1);
                    Z2(i,j)=B(2);
                end
     end
end

[X1,Y1,ZZ]=sph2cart(X,Y,Z1);
surf(X1,Y1,ZZ)
colormap('spring')

figure
[X1,Y1,ZZZ]=sph2cart(X,Y,Z2);
surf(X1,Y1,ZZZ)
fprintf('FINAAAAAL\n');

Answer 1

[SaiDileep Kola]

From the posted code I see Z1 and Z2 are always the same since B is always 1*1 double and only if condition(following B=real(sol(A))) is hit always, I also see that vpasolve is giving 4 different solutions, you may tweak to make use of all those solutions if they are the intended.