Since your problem is small, you're probably best off trying
root. Both of these are interfaces to MINPACK and call
HYBRJ. Since calculating a Jacobian matrix for your system shouldn't be hard (either do it by hand, or use your favorite computer algebra system, like SymPy, Sage, Maple, or Mathematica), you should supply a Jacobian matrix. The
HYBR functions in MINPACK use "Powell's hybrid method", which uses Newton's method and checks if Newton steps will be descent steps by comparing against least-squares minimization (specifically, does the Newton step also decrease the sum-of-squared-residuals). If Newton steps are not descent steps, then it falls back to the gradient of the sum-of-squared-residuals.
I think you're probably better off ignoring all of the alternatives for now. All of these other alternatives, broadly speaking, try to find low-dimensional structure in high-dimensional problems. Your problem has dimension 3, so the comparative savings of finding lower-dimensional structure aren't that great compared to, say, PDE discretizations that might have tens of thousands of variables.