import scipy.optimize as opt
from numpy import exp
def f(variables) :
(x,y) = variables
first_eq = x + y**2 - 4
second_eq = exp(x) + x*y - 3
return [first_eq, second_eq]
solution = opt.fsolve(f, (0.1, 1)) # fsolve(equations, X_0)
print(solution)
# [ 0.62034452 1.83838393]