How can I avoid a singular Jacobian when solving a system of equations with FindRoot?

I am trying to solve to following equation with findRoot.

rho0 = 1.0; rho1 = 1.0; rho2 = 10^(-10);
fr =
{-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x) == 0,
-(rho0 rho2)/(y^6) + (rho0 rho2)/y + (rho1 rho2)/(y - x) == 0};
frprecise = SetPrecision[fr, 32]
sol =
FindRoot[frprecise,{{x,0.015}, {y,0.8}},
AccuracyGoal -> 20,
PrecisionGoal -> 20,
WorkingPrecision -> 2 $MachinePrecision, MaxIterations -> 2000000] {-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /. sol No matter what values I take for the initial points, I still get a singular Jacobian. Would somebody kindly tell me how can I avoid the singular Jacobian and actually solve tis system to low residues (less than 10^(-30))? 2 Answers I believe that you cannot avoid a singular Jacobian because the Jacobian is singular at the solution you find. However, you can obtain more precision by using exact numbers for the inputs and by setting a higher working precision. rho0 = 1; rho1 = 1; rho2 = 10^(-10); fr = {-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x) == 0, -(rho0 rho2)/(y^6) + (rho0 rho2)/y + (rho1 rho2)/(y - x) == 0}; sol = FindRoot[fr, {{x, 15/1000}, {y, 8/10}}, WorkingPrecision -> 50, PrecisionGoal -> 30, MaxIterations -> 1000] {-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /. sol {-2.3492174793046040979*10^-30} Both NSolve[N@fr, {x, y}] and NSolve[frprecise, {x, y}, WorkingPrecision -> 32] give 30 complex solutions and no real solutions. To use FindRoot to search for a real solution is highly unlikely to be fruitful. The error (singular Jacobian) in FindRoot arises from FindRoot pursuing a zero by going out toward infinity. The magnitude of the functions decrease roughly as 1/Norm[{x, y}]. Eventually a point is reached at which numerical error diminishes the rank of the Jacobian from 2 to 1. sol = FindRoot[frprecise, {{x, 0.015}, {y, 0.8}}, AccuracyGoal -> 20, PrecisionGoal -> 20, WorkingPrecision -> 2$MachinePrecision,
MaxIterations -> 2000000]
(*  {x -> 1., y -> 591.241}  (FindRoot::jsing error omitted) *)

jac = D[frprecise /. Equal -> Subtract, {{x, y}}];
MatrixRank[jac /. sol]
(*  1  *)

Note the solution is reported in MachinePrecision even though higher precision was requested (maybe it's a bug). One can get the requested high precision by giving high precision initial values. Here I up the precision to 64, to get a residual less the the OP's requested goal.

frprecise = SetPrecision[fr, 64];
sol = FindRoot[frprecise, {{x, 0.01564}, {y, 0.864}},
AccuracyGoal -> 20, PrecisionGoal -> 20,
WorkingPrecision -> 4 \$MachinePrecision, MaxIterations -> 2000000]
(*  {x -> 0.9999999999999999..., y -> 3.5708094605541272691...*10^26}  (error omitted) *)

jac = D[frprecise /. Equal -> Subtract, {{x, y}}];
RowReduce[jac /. sol]
(*  {{1, 1.5685437702672624548...*10^-64}, {0, 0}}  *)

{-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /. sol
(*  {-2.80049*10^-37}  *)

We can get a much smaller residual by hand by going out further toward infinity, which does not mean it's getting closer to a root:

{-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /.
{x -> 1, y -> 10^100}
(*  {-1.*10^-110}  *)