Короткий ответ: используйте fsolve

, Как упомянуто в других ответах простое решение конкретной проблемы, которую Вы создали, должен использовать что-то как fsolve:

from scipy.optimize import fsolve
from math import exp

def f(vars):
    x, y = vars
    eq1 = x+y**2-4
    eq2 = exp(x) + x*y - 3
    return [eq1, eq2]

x, y =  fsolve(equations, (1, 1))

print(x, y)

Вывод:

0.6203445234801195 1.8383839306750887

Аналитические решения?

Вы говорите, как "решить", но существуют различные виды решения. Так как Вы упоминаете SymPy, я должен указать на самое большое различие между тем, что это могло средний, который является между аналитичен и числовой решения. Конкретным примером, который Вы дали, является тот, который не имеет (легкого) аналитического решения, но другие системы нелинейных уравнений делают. Когда существуют легко доступные аналитические решения, SymPY может часто находить их для Вас:

In [29]: from sympy import *                                                                                                                                  

In [30]: x, y = symbols('x, y')                                                                                                                               

In [31]: eq1 = Eq(x+y**2, 4)                                                                                                                                  

In [32]: eq2 = Eq(x**2 + y, 4)                                                                                                                                

In [33]: solve([eq1, eq2], [x, y])                                                                                                                            
Out[33]: 
⎡⎛ ⎛  5   √17⎞ ⎛3   √17⎞    √17   1⎞  ⎛ ⎛  5   √17⎞ ⎛3   √17⎞    1   √17⎞  ⎛ ⎛  3   √13⎞ ⎛√13   5⎞  1   √13⎞  ⎛ ⎛5   √13⎞ ⎛  √13   3⎞  1   √13⎞⎤
⎢⎜-⎜- ─ - ───⎟⋅⎜─ - ───⎟, - ─── - ─⎟, ⎜-⎜- ─ + ───⎟⋅⎜─ + ───⎟, - ─ + ───⎟, ⎜-⎜- ─ + ───⎟⋅⎜─── + ─⎟, ─ + ───⎟, ⎜-⎜─ - ───⎟⋅⎜- ─── - ─⎟, ─ - ───⎟⎥
⎣⎝ ⎝  2    2 ⎠ ⎝2    2 ⎠     2    2⎠  ⎝ ⎝  2    2 ⎠ ⎝2    2 ⎠    2    2 ⎠  ⎝ ⎝  2    2 ⎠ ⎝ 2    2⎠  2    2 ⎠  ⎝ ⎝2    2 ⎠ ⎝   2    2⎠  2    2 ⎠⎦

Примечание, что в этом примере SymPy находит все решения и не должен быть дан первоначальную оценку.

Точность числовых решений

Однако большинство систем нелинейных уравнений не будет иметь подходящего аналитического решения, настолько использующий SymPy как выше является большим, когда он работает, но не обычно применимый. Именно поэтому мы заканчиваем тем, что искали числовые решения даже при том, что с числовыми решениями: 1) У нас нет гарантии, что мы нашли все решения или "правильное" решение, когда существуют многие. 2) Мы должны обеспечить исходное предположение, которое не всегда легко.

признававший, что мы хотим числовые решения, что-то как fsolve будет обычно делать все, которое Вам нужно. Для этого вида проблемы SymPy, вероятно, будет намного медленнее, но он может предложить что-то еще, что находит (числовые) решения более точно:

In [50]: from sympy import *                                                                                                                                  

In [51]: x, y = symbols('x, y')                                                                                                                               

In [52]: nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1])                                                                                           
Out[52]: 
⎡0.620344523485226⎤
⎢                 ⎥
⎣1.83838393066159 ⎦

In [53]: nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1], prec=50)                                                                                  
Out[53]: 
⎡0.62034452348522585617392716579154399314071550594401⎤
⎢                                                    ⎥
⎣ 1.838383930661594459049793153371142549403114879699 ⎦