sistema de dos ecuaciones que para resolverlo utilizo un método numérico

from scipy.optimize import fsolve

def equations(p):
    x, y = p
    return (y - x**2 + 1, np.exp(x) + x - x*y - 1)

x1, y1 =  fsolve(equations, (1, 1))

print(x1, y1)