Finding the roots of a system of non-linear equations that has multiple roots with python fsolve

I have the following system of non-linear equations which I want to find its roots:

x - exp(a x + b y) = 0
y - exp(c x + d y) = 0 The code I'm using to find its roots is:

equations = lambda x, kernel: np.array([x[0] - np.exp(kernel[0] * x[0] + kernel[2] * x[1]), x[1] - np.exp(kernel[1] * x[1] + kernel[3] * x[0])])

kernels = np.array([kernel0, kernel1, kernel2, kernel3])
x_init = np.array([x_init0, x_init1])
x_sol = fsolve(two_equations, x_init, args=(kernels))

From the equations I know that this system in some situations has two answers for each variable: (x_sol1, x_sol2) and (y_sol1, y_sol2).

Is there a clean way to pass multiple initial guesses to this fsolve function to get both roots for each variable? (instead of using a for loop) I know how to do it for a system of one equation only but I couldn't use that method for this case.

Solution

You can reduce the system to a single univariate equation by eliminating y.

y = (ln(x) - a x) / b

so that

(ln(x) - a x) / b - exp(c x + d (ln(x) - a x) / b) = 0