# Solve two equation by fixed point iteration method

I am tryin to reproduce the figure.2 of this paper. In order to obtain this plot, I need to find the global minimum of equation.(7) or equivalently solving two coupled equations of $$\frac{\partial \Omega}{\partial M}=0$$ and $$\frac{\partial \Omega}{\partial q}=0$$. Finding the global minimum of the equation.(7) is taking much time, specially when I consider other parameters as variables. So I prefer to find the solution of two equations $$\frac{\partial \Omega}{\partial M}=0$$ and $$\frac{\partial \Omega}{\partial q}=0$$ using the fixed-point iteration method. To this purpose, at each level of iteration, I solve the equation $$\frac{\partial \Omega}{\partial M}=0$$ and find the minimum of $$\Omega$$ with respect to the $$q$$. Here, is the definition of $$\Omega$$

Omega[m_?NumberQ,m0_?NumberQ,b_?NumberQ,Lambda_?NumberQ,Mu_?NumberQ,q_?NumberQ]:=(m - m0)^2*(Lambda^2)/(4*b)
- 3/(Pi^2) * Quiet[NIntegrate[ p*((Mu - Sqrt[p^2 + (Sqrt[k^2 + m^2] +q/2)^2])*
HeavisideTheta[Mu - Sqrt[p^2 + (Sqrt[k^2 + m^2] + q/2)^2]] + (Mu - Sqrt[p^2 + (Sqrt[k^2 + m^2] - q/2)^2])*
HeavisideTheta[Mu - Sqrt[p^2 + (Sqrt[k^2 + m^2] -q/2)^2]]),{p, 0.0000001, 100}, {k, 0.0000001, 100},
MinRecursion -> 5, PrecisionGoal -> 12,AccuracyGoal -> 5]]+ 3/(4*Pi^(5/2)) Quiet[NIntegrate[ 1/t^(5/2)*
(Exp[-(Sqrt[k^2 + m^2] + q/2)^2*t] + Exp[-(Sqrt[k^2 + m^2] - q/2)^2*t]), {k, 0.0000001,100},
{t,Lambda^(-2), 100}, MinRecursion -> 5, PrecisionGoal -> 12, AccuracyGoal -> 5]];


I arranged the $$\frac{\partial \Omega}{\partial m}=0$$ to the form $$m=GM(m,m0,b,Lambda,Mu,q)$$, which is,:

GM[m_?NumberQ,m0_?NumberQ,b_?NumberQ,Lambda_?NumberQ,Mu_?NumberQ,q_?NumberQ]:= -(2*b/Lambda^2)*
(3*m/Pi^2 * Quiet[NIntegrate[ p/Sqrt[k^2+m^2]*((Sqrt[k^2 + m^2] + q/2)/Sqrt[p^2+(Sqrt[k^2 + m^2]+q/2)^2]
*HeavisideTheta[Mu - Sqrt[p^2 + (Sqrt[k^2 + m^2] + q/2)^2]]
+ (Sqrt[k^2 + m^2] - q/2)/Sqrt[p^2 + (Sqrt[k^2 + m^2] - q/2)^2]
*HeavisideTheta[Mu - Sqrt[p^2 + (Sqrt[k^2 + m^2] - q/2)^2]]),{p, 0.0000001, 100},{k,0.0000001,100},MinRecursion ->5, PrecisionGoal->12,AccuracyGoal->5]]
-(3*m)/(2 Sqrt[Pi^5]) Quiet[NIntegrate[1/t^(3/2)*((Sqrt[k^2 + m^2] + q/2)/Sqrt[k^2 + m^2]*Exp[-(Sqrt[k^2 + m^2] + q/2)^2*t] + (Sqrt[k^2 + m^2] - q/2)/Sqrt[k^2 + m^2]*Exp[-(Sqrt[k^2 + m^2] - q/2)^2*t]), {k, 0.0000001,100}, {t,Lambda^(-2),100}, MinRecursion -> 5, PrecisionGoal->12,AccuracyGoal->5]])


I defined a function to find the minimum of $$\Omega$$ at each iteration for the obtained value of $$m$$:

func[m_?NumberQ, m0_?NumberQ,b_?NumberQ,Lambda]_?NumberQ,Mu_?NumberQ] := Module[{startpoint=
Omega[m,m0,b,Lambda,Mu,0]},Monitor[Do[ If[Omega[m, m0, b,Lambda,Mu,q]<=startpoint,qu = q;startpoint=Omega[m,m0,b,Lambda,Mu,q]];,{q, 0, 0.8, 0.001}], q];Return[qu];]


With all these the main code is:

IniM = 0.4; IniQ = 0.4;DenErr = 10^(-20);Mu = 0;
Block[{Lambda = 0.86, b = 6, m0 = 0},
Monitor[ Do[ error = 1; SecM = 0; SecQ = 0; LoopCount = 0;
While[error > 0.001, LoopCount += 1; SecM = GM[IniM, m0, b, Lambda, Mu, IniQ];
SecQ = func[SecM, m0, b, Lambda, Mu];
error = Sqrt[(SecM - IniM)^2/(Max[SecM, IniM])^2 + (SecQ-IniQ)^2/((Max[SecQ, IniQ])^2 + DenErr)];
IniM = SecM; IniQ = SecQ; Print[IniM, " ", IniQ, " ", error, " ", LoopCount];];
Print[Mu, " ", IniM, " ", IniQ, " ", LoopCount]; Mu += 0.005;,{i, 0, 100, 1}], i]]


The question is, with the proper choice of initial value of iteration, I am able to obtain the result up to the $$\mu = 0.415$$. However, for the loop $$i=84$$, or $$\mu=0.42$$, the iteration fail to give an answer after hours of running. I would be grateful if someone could help me to modify the code in a way this problem can be solved and get results for $$i \geq 84$$ or $$\mu \geq0.42$$.

P.S: at the low value of $$\mu$$ the solution for $$q$$ is zero. so I added a small value to the denominator of error, for the part related to $$q$$ to avoid divergence, to avoid divergence.