Rootfinding (system of nonlinear equations) with bound constraints

Question

I would like to solve in Matlab a system of nonlinear equations with bound constraints, i.e. lb<=x<=ub. Matlab suggests to convert the problem to a minimization and use lsqnonlin or fmincon. I want instead to use a transformation with fsolve. As an example, I picked this one: https://www.mathworks.com/help/optim/ug/nonlinear-systems-with-constraints.html

This is a system of 2 nonlinear equations in two unknowns, with 4 different solutions. I impose the bounds x>=0 and only one solution satisfies the bounds, namely x*=(10,20). To solve the system

F(x)=0 s.t. x>=0,

I call fsolve with F(z.^2) and I find the solution z* and transform it back to x=z.^2. The square transformation makes sure that the variable x never becomes negative. As a starting point I choose a feasible value, avoiding x0=(0,0).

However, the transformation method does not work. I also tried lsqnonlin with bounds and it does not work as well. It works only if I choose an initial condition very close to the true solution (10,20), but suppose in a more complicated problem I don't know what the solution is. Q: Should I use other transformations, such as x=exp(z)? Does the transformation have to be a one-to-one function (the square is obviously not)? Can someone help me on this? Thanks a lot!

MINIMUM WORKING EXAMPLE HERE:

%% fsolve with bound constraints
% Reference: https://www.mathworks.com/help/optim/ug/nonlinear-systems-with-constraints.html
% The system of equations have 4 different solutions but only one obeys the
% nonnegativity constraint x(1)>=0, x(2)>=0
% (-1,-2) NOT FEASIBLE
% (10,-2) NOT FEASIBLE
% (-1,20) NOT FEASIBLE
% (10,20) NOT FEASIBLE
% Transformation: instead of solving F(x)=0 where x is unbounded, solve for
% F(z.^2)=0 where z is unbounded but of course z.^2>=0. 
clear
clc
close 
x0 = [5,10];
%opts = optimoptions(@fsolve,'Display','off');
%% Transformation method
x0 = sqrt(x0);
[x_star,~,flag] = fsolve(@fbnd_trans,x0); 
x_star = x_star.^2;
if flag<=0
    warning('Equations not solved!')
end
disp('Solution x* = ')
disp(x_star)
disp('Residuals = ')
Residuals = 
disp(fbnd(x_star))
 
%% lsqnonlin
lb = [0,0]';
ub = [inf,inf]';
[x_star,~,~,flag] = lsqnonlin(@fbnd,x0,lb,ub);
if flag<=0
    warning('Equations not solved!')
end
disp('Solution x* = ')
disp(x_star)
disp('Residuals = ')
disp(fbnd(x_star))

%------------------------- SUBFUNCTIONS ----------------------------------% 
function F = fbnd(x)
F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end
function F = fbnd_trans(x)
% Impose x>=0 bound
x = x.^2;
F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end