solving a nonlinear set of equations while have a dependency in it

I have the following set of 7 equations. There is a dependency between first 6 equations as following:

(Z3^2/Z1)-Z2=(X1^2/X3)-X2

and so I think T have, in fact, 5 equations and that's why I added the last equation as well (to have 6 unknowns and 6 independent eqs.) I solved this system with Mathematica and I know that it should return 2 answers as following:

{a1 -> 2.99957, a2 -> 1.99, a3 -> 0.999892, b1 -> 0.178447, b2 -> -0.19892, b3 -> 0.124498},
{a1 -> 2.99957, a2 -> 2.00, a3 -> 0.999892, b1 -> 0.199952, b2 -> 0.198922, b3 -> 0.199778}

However, in Matlab I couldn't get a solution using "solve". I used the following code; would you please have a look on it and let me know what is the problem (Should I define for "solve" that there is a dependency between first 6 eqs., etc)

Since based on constraint I have in my problem, the solution method must be numerical (and I should be able to adjust the error band) so, "fsolve" is also good here if can return both of the answers for me. (Thanks a lot)

syms a1 a2 a3 b1 b2 b3

X1=0.1866667;
X2=1.9866667;
X3=0.9866667;
Z1=2.96;
Z2=1.96;
Z3=0.16;

E=[a1 b3 b2;                            % all are my 6 unknowns %
   b3 a2 b1;
   b2 b1 a3];
L=eig(E);
s=eval(4/3*pi/sqrt(L(1)*L(2)*L(3)));    % a volume equation based on unknowns %    

[a1, a2, a3, b1, b2, b3]=solve(b1-b2*b3/a1==X1,...
                               a2-b3^2/a1==X2,...   
                               a3-b2^2/a1==X3,...   
                               a1-b2^2/a3==Z1,...   
                               a2-b1^2/a3==Z2,...   
                               b3-b1*b2/a3==Z3,...
                               s==1.743);                           
a1=eval(a1)
a2=eval(a2)
a3=eval(a3)
b1=eval(b1)
b2=eval(b2)
b3=eval(b3)