Solving mean field nonlinear equations uisng Findroot

I am trying to solve three mean field equations using Findroot.
I tried several options but no success.
Could someone help me out and code is attached/
Thank you

meanfieldeqn.nb (55.8 KB)

The functions that you are integrating are continuous, but just barely. They contain the expression

Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] - Sqrt[
     a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])] + 
 Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] + Sqrt[
     a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])]

that is almost the same as

2 UnitStep[\[Epsilon] - (a^2 + 4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)]

You may use this fact to simplify the functions.

Thank you for the reply.

The one you gave me was for sum of Tanh functions.

Could you please give me an expresssion for difference of Tanh functions as well?

Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] - Sqrt[
     a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])]- 
 Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] + Sqrt[
     a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])]

I would say

2 UnitStep[\[Epsilon] - (a^2 + 
       4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)] - 2

To check it:

Manipulate[
 Plot[Tanh[
    2500 (\[Epsilon] + \[Epsilon]f - \[Mu] - 
       Sqrt[a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])] - 
   Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] + 
       Sqrt[a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])] - (
     2 UnitStep[\[Epsilon] - (a^2 + 
           4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)] - 
     2), {\[Epsilon], -1, 1}],
 {{\[Mu], -.2152}, -1, 1},
 {{\[Epsilon]f, 1.2474}, 0, 2},
 {{a, .1}, -1, 1}]

Thank you again.
I replaced sum and diff tanh functions with unitstep function but still solutions didn’t satisfy nonlinear equations.

In[1]:= d = 1;
T = 0.0001 d;
\[Epsilon]f0 = - 0.102001 d;
p = 0.875;
q0 = 0.5;
v = 0.5 d;
ds[\[Epsilon]_] := 3/(4 Sqrt[2] d) Sqrt[(\[Epsilon] + d)/d];

(*we calculate renormalized position of the f-level, chemical \
potential and slave boson amplitude by solving three nonlinear \
equations*)

\[Zeta][\[Epsilon]_, \[Mu]_] := \[Epsilon] - \[Mu];  \
(*\[Epsilon]=k^2/(2 m)-d*)
Rk[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
 Sqrt[(\[Epsilon]f - \[Zeta][\[Epsilon], \[Mu]])^2 + (2 v a)^2];
Enkp[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (\[Epsilon]f + \[Zeta][\[Epsilon], \[Mu]] + 
     Rk[\[Epsilon], \[Mu], \[Epsilon]f, a]);
Enkm[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (\[Epsilon]f + \[Zeta][\[Epsilon], \[Mu]] - 
     Rk[\[Epsilon], \[Mu], \[Epsilon]f, a]);

In[11]:= fwp[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (1 - Tanh[Enkp[\[Epsilon], \[Mu], \[Epsilon]f, a]/(2 T)]);
fwm[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (1 - Tanh[Enkm[\[Epsilon], \[Mu], \[Epsilon]f, a]/(2 T)]);
diff[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1 - UnitStep[\[Epsilon] - (a^2 + 
       4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)];
sum[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1 - UnitStep[\[Epsilon] - (a^2 + 
       4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)];
eq1int[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  diff[\[Epsilon], \[Mu], \[Epsilon]f, a]/
  Rk[\[Epsilon], \[Mu], \[Epsilon]f, a];
eq2int[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, 
   a_] := (\[Zeta][\[Epsilon], \[Mu]] - \[Epsilon]f) \
eq1int[\[Epsilon], \[Mu], \[Epsilon]f, a];

eq3int[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  sum[\[Epsilon], \[Mu], \[Epsilon]f, a];

fint1[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  ds[\[Epsilon]] eq1int[\[Epsilon], \[Mu], \[Epsilon]f, a];
fint2[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  ds[\[Epsilon]] eq2int[\[Epsilon], \[Mu], \[Epsilon]f, a];
fint3[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  ds[\[Epsilon]] eq3int[\[Epsilon], \[Mu], \[Epsilon]f, a];

In[21]:= eqn1[\[Mu]_?NumericQ, \[Epsilon]f_?NumericQ, 
  a_?NumericQ] := \[Epsilon]f - \[Epsilon]f0 - 
  v^2 NIntegrate[
    fint1[\[Epsilon], \[Mu], \[Epsilon]f, a], {\[Epsilon], -d, d}, 
    MaxRecursion -> 20, 
    Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]
eqn2[\[Mu]_?NumericQ, \[Epsilon]f_?NumericQ, a_?NumericQ] := 
 2 (q0 - a^2 ) - p - 
  NIntegrate[
   fint2[\[Epsilon], \[Mu], \[Epsilon]f, a], {\[Epsilon], -d, d}, 
   MaxRecursion -> 20, 
   Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]
eqn3[\[Mu]_?NumericQ, \[Epsilon]f_?NumericQ, a_?NumericQ] := 
 p - NIntegrate[
   fint3[\[Epsilon], \[Mu], \[Epsilon]f, a], {\[Epsilon], -d, d}, 
   MaxRecursion -> 20, 
   Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]

In[24]:= mysol = 
 FindRoot[{eqn1[\[Mu], \[Epsilon]f, a], eqn2[\[Mu], \[Epsilon]f, a], 
   eqn3[\[Mu], \[Epsilon]f, a]}, {{\[Mu], -0.42628}, {\[Epsilon]f, 
    0.035940}, {a, -0.37342}}, WorkingPrecision -> 16, 
  AccuracyGoal -> 15, MaxIterations -> Infinity]


During evaluation of In[24]:= FindRoot::jsing: Encountered a singular Jacobian at the point {\[Mu],\[Epsilon]f,a} = {0.845381,0.196724,-0.730453}. Try perturbing the initial point(s).

Out[24]= {\[Mu] -> 0.845381, \[Epsilon]f -> 0.196724, a -> -0.730453}

In[25]:= {eqn1[\[Mu], \[Epsilon]f, a], eqn2[\[Mu], \[Epsilon]f, a], 
  eqn3[\[Mu], \[Epsilon]f, a]} /. mysol

Out[25]= {0.0630848, -0.292974, -0.125}

The formulas for the sum and difference of tanh are not the same. Also, they are only approximations to be estimated.