---

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, root

f1 = lambda x1: 1/(x1+1.5) - 1
f2 = lambda x1: (1.2*x1)**3 - 1

x = np.linspace(-1, 2, 500)
plt.plot(x, f1(x), label='x2 = 1/(x1+1.5)-1')
plt.plot(x, f2(x), label='x2 = (1.2*x1)^3-1')
plt.xlabel('x1'); plt.ylabel('x2'); plt.grid(); plt.legend()

def sys(v):
    return [v[1] - f1(v[0]), v[1] - f2(v[0])]

guess = [0.0, 0.0]
sol_fsolve = fsolve(sys, guess)
sol_root = root(sys, guess).x

print(f"fsolve: x1={sol_fsolve[0]:.6f}, x2={sol_fsolve[1]:.6f}")
print(f"root:   x1={sol_root[0]:.6f}, x2={sol_root[1]:.6f}")

plt.plot(sol_fsolve[0], sol_fsolve[1], 'ro', markersize=8, label='Решение')
plt.legend(); plt.show()