Решение системы с помощью Sympy:
from sympy import symbols, Eq, solve
x, y = symbols("x y")
equations = [
Eq(x**2 + 5, y),
Eq(y**2 + x, 3)
]
results = solve(equations)
print(results)
# Вывод: [{x: 3 - (sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2 - sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2)**2, y: sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2 - sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2}, {x: 3 - (sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2 + sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2)**2, y: sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2 + sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2}]
# Вывод результатов в более читаемом виде
for result in results:
print('x =', result[x])
print(f'({result[x].evalf()})')
print()
print('y =', result[y])
print(f'({result[y].evalf()})')
print('\n')
# Вывод:
# x = 3 - (sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2 - sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2)**2
# (-0.136141938982763 - 1.78387663318005*I)
# y = sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2 - sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2
# (1.83631878514418 + 0.485720847494351*I)
# x = 3 - (sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2 + sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2)**2
# (-0.136141938982763 + 1.78387663318005*I)
# y = sqrt(8 - 2*(209/16 + sqrt(235023)*I/144)**(1/3) + 2/sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3)) - 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)))/2 + sqrt(4 + 34/(3*(209/16 + sqrt(235023)*I/144)**(1/3)) + 2*(209/16 + sqrt(235023)*I/144)**(1/3))/2
# (1.83631878514418 - 0.485720847494351*I)
Видим, что система уравнений из вопроса имеет только комплексные решения. Правда, WolframAlfa дает для этой системы на два решения больше (кликабельно):
Если комплексные решения не нужны, проверяйте свойство is_real
(например result[x].evalf().is_real
), и решения где это свойство равно False не учитывайте:
# Убрал +5 в первом уравнении, чтобы появились действительные корни:
equations = [
Eq(x**2, y),
Eq(y**2 + x, 3)
]
...
for result in results:
if result[x].evalf().is_real and result[y].evalf().is_real:
print('x =', result[x])
print(f'({result[x].evalf()})')
print()
print('y =', result[y])
print(f'({result[y].evalf()})')
print('\n')
Вывод:
x = 3 - (-sqrt(-2*(sqrt(257)/16 + 129/16)**(1/3) - 8/(sqrt(257)/16 + 129/16)**(1/3) + 2/sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3)) + 8)/2 + sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3))/2)**2
(1.16403514028977)
y = -sqrt(-2*(sqrt(257)/16 + 129/16)**(1/3) - 8/(sqrt(257)/16 + 129/16)**(1/3) + 2/sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3)) + 8)/2 + sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3))/2
(1.35497780782942)
x = 3 - (sqrt(-2*(sqrt(257)/16 + 129/16)**(1/3) - 8/(sqrt(257)/16 + 129/16)**(1/3) + 2/sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3)) + 8)/2 + sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3))/2)**2
(-1.45262687883384)
y = sqrt(-2*(sqrt(257)/16 + 129/16)**(1/3) - 8/(sqrt(257)/16 + 129/16)**(1/3) + 2/sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3)) + 8)/2 + sqrt(8/(sqrt(257)/16 + 129/16)**(1/3) + 4 + 2*(sqrt(257)/16 + 129/16)**(1/3))/2
(2.11012484911056)