Problem with FindRoot function

Question

I have set up a code to solve simultaneous integral equations. These result essentially from taking integrals of the Fermi-Dirac distribution. Although there are no obvious problems with the code it does not give the correct answers for the roots. I have checked that in the limiting case when I could evaluate the integrals analytically and find the roots exactly. Here is the code :

dv1 = 1
dv2 = 1
u1 = 5
u2 = 8
η1 = 1
η2 = 0.5
N1 = 2 * η1 * dv1
N2 = 2 * η2 * dv2
K = η1 * (dv1 ^ 3 / 12 + u1 ^ 2 * dv1) + η2 * (dv2 ^ 3 / 12 + u2 ^ 2 * dv2)
f[b_, x_, y_] := 1/Sqrt[2*η1]*NIntegrate[b^(-1/2)*z^(-1/2)*η1*Exp[-z + b*x*η1]/(1 + Exp[-z + b*x*η1] + Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}]
g[b_, x_, y_] := 1/Sqrt[2*η2]*NIntegrate[b^(-1/2)*z^(-1/2)*η2*Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] + Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]
h[b_, x_, y_] := 
 2^(-1/2)*(η1)^(-3/2)*NIntegrate[b^(-3/2)*z^(1/2)*η1*Exp[-z + b*x*\[Eta]1]/(1 + Exp[-z + b*x*η1] + Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}] + 
  2^(-1/2)*(η2)^(-3/2)*NIntegrate[b^(-3/2)*z^(1/2)*η2*Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] + Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]
sol = FindRoot[{N1 - f[b, x, y] == 0, N2 - g[b, x, y] == 0, K - h[b, x, y] == 0}, {{b, 0.0349599}, {x, -54.4214}, {y, -128.67}}, WorkingPrecision -> 40] /. Rule -> List
b = Re[sol[[1]][[2]]]
x = Re[sol[[2]][[2]]]
y = Re[sol[[3]][[2]]]

I am trying to solve the equations for the variables b,x,y. The problem seems to be with the FindRoot part of the code. Even if I specify as the initial values for these variables the values that are the true solutions I do not get them in return. As I am only a beginner in Mathematica I don't have an idea how to proceed now. There might be more than one root to the system of equations but I don't know for sure.

Answer 1

Need to change N1 to n1 as N is a built-in function name. Also need to change function arguments to x_?NumericQ as I did below. Then change 0.5 to 1/2 to set up arbitrary-precision calculations. Then this seems to work:

dv1 = 1
dv2 = 1
u1 = 5
u2 = 8
η1 = 1
η2 = 1/2
n1 = 2*η1*dv1
n2 = 2*η2*dv2
K = η1*(dv1^3/12 + u1^2*dv1) + η2*(dv2^3/12 + u2^2*dv2)
f[b_?NumericQ, x_?NumericQ, y_?NumericQ] := 
 1/Sqrt[2*η1]*
  NIntegrate[
   b^(-1/2)*z^(-1/2)*η1*
    Exp[-z + b*x*η1]/(1 + Exp[-z + b*x*η1] + 
       Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}]

g[b_?NumericQ, x_?NumericQ, y_?NumericQ] := 
 1/Sqrt[2*η2]*
  NIntegrate[
   b^(-1/2)*z^(-1/2)*η2*
    Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] + 
       Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]

h[b_?NumericQ, x_?NumericQ, y_?NumericQ] := 
 2^(-1/2)*(η1)^(-3/2)*
   NIntegrate[
    b^(-3/2)*z^(1/2)*η1*
     Exp[-z + b*x*η1]/(1 + Exp[-z + b*x*η1] + 
        Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}] + 
  2^(-1/2)*(η2)^(-3/2)*
   NIntegrate[
    b^(-3/2)*z^(1/2)*η2*
     Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] + 
        Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]

sol = FindRoot[{n1 - f[b, x, y] == 0, n2 - g[b, x, y] == 0, 
    K - h[b, x, y] == 0}, {{b, 
     0.0349599}, {x, -54.4214}, {y, -128.67}}, 
   WorkingPrecision -> 40] /. Rule -> List
b = Re[sol[[1]][[2]]]
x = Re[sol[[2]][[2]]]
y = Re[sol[[3]][[2]]]
theSol = #[[2]] & /@ sol
Apply[{(n1 - f[#1, #2, #3]), n2 - g[#1, #2, #3], 
   K - h[#1, #2, #3]} &, theSol]*emphasized text*

{0.,3.33067*10^-16,1.42109*10^-14}

Answer 2

We have the three simultaneous equations:

$$ \begin{align} n1-f(b,x,y)&=0\\ n2-g(b,x,y)&=0\\ K-h(b,x,y)&=0 \end{align} $$

Now when I run the code above, I get:

{0.04331798263222832287097452988786750625601,
-14.02830453956710576253245993437615631186,
-47.91130641078542149075012315062585520063}

and back-substituting those values into the three equations returns:

{0., 3.33067*10^-16, 1.42109*10^-14} 

However, when I back-substitute your solutions:

{0.0349599,-54.4214,-128.67} 

I do not get answers close to zero:

{0.841642, 0.432857, 26.3829}

so must conclude the problem is not with FindRoot but rather before it, that is, the integrations perhaps but do not receive any type of error message with the integrations.