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.