Suppose I have a multivariate function, f(x, y), with x \in (0, 1) and y \in (0, 1). I want to find the roots of following equation

f(x, y) - c = 0

where c = [c_1\ c_2]^T is a constant vector. In my case, there are `log`

functions in f.

I have tried `NLsolve`

. It can find correct root sometimes, but it strongly depends on how close the initial guess to the root. If a bad initial guess is given, the `log`

function will complain for receiving a negative value. I wonder if I can provide the bounds for x and y to `NLsolve`

?

I have also tried `NonlinearSolve`

with default settings. It always fails.

Can someone provide me some guide on how to find roots of such system? BTW, the analytical expression of f may not know and sometimes extremely costly to compute.

MWE below

```
using StaticArrays
# a typical set of parameters
param = (1.0, 0.2, 0.01, 10.0, 40.0, 40.0)
# the objective function
f = (x, y) -> μ(x, y, param)
# an example of c can be obtained from below
c = f(0.2, 0.4)
# try to find the roots of
g = (x, y) -> f(x, y) - c # should give at least one root which is [0.2, 0.4]. Other roots are possible.
function γA(ϕA, ϕB, param)
αA, αB, αS, χABN, χASN, χBSN = param
ϕS = 1 - ϕA - ϕB
return (1+log(ϕA))/αA + χABN*ϕB + χASN*ϕS
end
function γB(ϕA, ϕB, param)
αA, αB, αS, χABN, χASN, χBSN = param
ϕS = 1 - ϕA - ϕB
return (1+log(ϕB))/αB + χABN*ϕA + χBSN*ϕS
end
function γS(ϕA, ϕB, param)
αA, αB, αS, χABN, χASN, χBSN = param
ϕS = 1 - ϕA - ϕB
return (1+log(ϕS))/αS + χASN*ϕA + χBSN*ϕB
end
function μA(ϕA, ϕB, param)
return γA(ϕA, ϕB, param) - γS(ϕA, ϕB, param)
end
function μB(ϕA, ϕB, param)
return γB(ϕA, ϕB, param) - γS(ϕA, ϕB, param)
end
function μ(ϕA, ϕB, param)
return SVector(μA(ϕA, ϕB, param), μB(ϕA, ϕB, param))
end
```