Cómo resolver un par de ecuaciones no lineales usando Python?

Si prefiere sympy, puede usar nsolve.

>>> nsolve([x+y**2-4, exp(x)+x*y-3], [x, y], [1, 1])
[0.620344523485226]
[1.83838393066159]

El primer argumento es una lista de ecuaciones, el segundo es una lista de variables y el tercero es una suposición inicial.

Prueba este, te aseguro que funcionará perfectamente.

    import scipy.optimize as opt
    from numpy import exp
    import timeit

    st1 = timeit.default_timer()

    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) )
    print(solution)


    st2 = timeit.default_timer()
    print("RUN TIME : {0}".format(st2-st1))

->

[ 0.62034452  1.83838393]
RUN TIME : 0.0009331008900937708

FYI. como se mencionó anteriormente, también puede usar 'Aproximación de Broyden 'reemplazando' fsolve 'por'broyden1'. Funciona. Lo hice.

No se exactamente como funciona la aproximación de Broyden, pero tomó 0.02 s.

Y le recomiendo que no utilice las funciones de Sympy

Puede usar el paquete openopt y su método NLP. Tiene muchos algoritmos de programación dinámica para resolver ecuaciones algebraicas no lineales que consisten en:
GoldenSection, scipy_fminbound, scipy_bfgs, scipy_cg, scipy_ncg, amsg2p, scipy_lbfgsb, scipy_tnc, bobyqa, ralg, ipopt, scipy_slsqp, scipy_cobyla, lincher, algencan, que puede elegir.
Algunos de los últimos algoritmos pueden resolver problemas de programación no lineal restringidos. Por lo tanto, usted puede introducir su sistema de ecuaciones a openopt.NLP () con una función como esta:

lambda x: x[0] + x[1]**2 - 4, np.exp(x[0]) + x[0]*x[1]

Conseguí que el método de Broyden funcionara para ecuaciones no lineales acopladas (generalmente que involucran polinomios y exponenciales) en IDL, pero no lo he probado en Python:

Http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.broyden1.html#scipy.optimize.broyden1

Scipy.optimizar.broyden1

scipy.optimize.broyden1(F, xin, iter=None, alpha=None, reduction_method='restart', max_rank=None, verbose=False, maxiter=None, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None, tol_norm=None, line_search='armijo', callback=None, **kw)[source]

Encuentra una raíz de una función, usando la primera aproximación jacobiana de Broyden.

Este método también se conoce como " Broyden's good método".

from scipy.optimize import fsolve

def double_solve(f1,f2,x0,y0):
    func = lambda x: [f1(x[0], x[1]), f2(x[0], x[1])]
    return fsolve(func,[x0,y0])

def n_solve(functions,variables):
    func = lambda x: [ f(*x) for f in functions]
    return fsolve(func, variables)

f1 = lambda x,y : x**2+y**2-1
f2 = lambda x,y : x-y

res = double_solve(f1,f2,1,0)
res = n_solve([f1,f2],[1.0,0.0])

Una alternativa a fsolve es root:

import numpy as np
from scipy.optimize import root    

def your_funcs(X):

    x, y = X
    # all RHS have to be 0
    f = [x + y**2 - 4,
         np.exp(x) + x * y - 3]

    return f

sol = root(your_funcs, [1.0, 1.0])
print(sol.x)

Esto imprimirá

[0.62034452 1.83838393]

Si luego marca

print(your_funcs(sol.x))

Obtienes

[4.4508396968012676e-11, -1.0512035686360832e-11]

, Confirmando que la solución es correcta.