I got Broyden's method to work for coupled non-linear equations (generally involving polynomials and exponentials) in IDL, but I haven't tried it in Python:
scipy.optimize.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]
Find a root of a function, using Broyden’s first Jacobian approximation.
This method is also known as “Broyden’s good method”.