for loop solve implicit equation with 2 roots

Hello, I want to find all the possible roots for the equation gammameta for given ranges of values for theta and r. You can see my trying in the script below. I have two problems with my script:

1. Each array in xsol is a 1x1 sym, but I know that for some cases I can obtain 2 roots. I cannot get why the programm is only showing me one.
2. When I run the vpa function , I do not obtain a numerical result. I still obtain sym

I don't know what I am doing wrong.

s=0.4;
rho=0.01;
b=2.5;
tau=0.286;

%ranges
theta_pool = 0.01:0.01:0.75;
r_pool = 0.02:0.001:0.035;
beta_pool= 0:0.001:0.2;
xsol=cell(75,16)
%Call fsolve in a loop
for i = 1 : numel(theta_pool)
    theta = theta_pool(i);
for j = 1: numel(r_pool)
       r = r_pool(j);
  syms x 
  gamameta= ((r*((1-tau)*s*(x^(1-s))*(1-theta))/(r-theta*(1-tau)*s*x^(1-s)))- rho -(x^(s)/b)-r+tau/b);
  xsol{i,j}=solve(gamameta==0,x);
end
end 
xsol2=cellfun(@(x)vpa(x),xsol,'un',false);