Using fzero to solve an equation with two different constants every time

Hello, I want to solve for x(i) the function F for multiple values of theta and r and record each x(i). I used the code below to do it and I obtain that all x(i) are copmplex numbers. However, when I specify values for theta (0.06) and r=0.026 and use fzero, I obtain x= 0.995. Note that the values of theta and r that I chose are inside the ranges I chose for those variables. I think I am doing the nested for loop wrong.

First, I show you the code with the nested loop and then the code for the specific values of theta and r.

Thank you for your help!

%NESTED LOOP CODE
function xsol=main() 
s=0.4;
rho=0.01;
b=2.5;
tau=0.286;
x0=[0.099;0.1];

%values for tau
theta_pool = 0.01:0.01:0.75;
r_pool = 0.02:0.0001:0.035;
%preallocate result matrix 
%x_all = zeros(numel(x0),numel(tau1_pool))

% switch display off
%options = optimoptions('fzero','Display','off');

%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);
    xsol(i,j) = fzero(@obj_fun,x0);
   end
end

function F = obj_fun(x)
F=(r*((1-tau)*s*(x^(1-s))*(1-theta))/(r-theta*(1-tau)*s*x^(1-s)))- rho -(x^(s)/b)-r+tau/b;
end
end
%call the result by typing result=main 


% fzero FOR SPECIFIC THETA AND R VALUES CODE

myfun=@(x,r,tau,s,theta,rho,b) (r*((1-tau)*s*(x^(1-s))*(1-theta))/(r-theta*(1-tau)*s*x^(1-s)))- rho -(x^(s)/b)-r+(tau/b);
r=0.026
tau=0.286
s=0.4
rho=0.01
theta=0.06
b=2.5
x0=[0.099; 0.1];
fun=@(x) myfun(x,r,tau,s,theta,rho,b);
xsol= fzero(fun,x0)