How does sympy's nsolve do what it does?
As @hpaulj mentioned in the comment, SymPy's nsolve uses mpmath.findroot. This is what happens (roughly speaking) when you execute nsolve(equations, variables, initial_guess) (I'm using SymPy 1.13.3 at the time of writing this answer):
the symbolic equation (or system of equations) are converted to a numerical function using lambdify.
If a system of equations is provided, the Jacobian is also computed. Then, it is converted to a numerical function using lambdify.
mpmath.findroot is executed, like this:
x = sympify(findroot(f, x0, **kwargs)) for a single equation. If no keyword arguments are provided to nsolve, then the secant method will be used.x = findroot(f, x0, J=J, **kwargs) for a system of equations. If no keyword arguments are provided to nsolve, then the multidimensional Newton method is used.Note that we can actually control the behavior of mpmath.findroot by providing appropriate keyword arguments to nsolve, which are then passed down to mpmath.findroot.
The result x is parsed in order to return a solution with the requested form (dict or list).
Note that we can force nsolve to use the multidimensional Newton method also for a single equation by wrapping the arguments in lists, like this: nsolve([single_equation], [variable], [initial_guess]). This trick can be used when the secant method fails to find a solution.