I have two equations:
$$ \frac{2}{3}\thinspace\left(\frac{\varepsilon}{\Delta}\right)^{\frac{3}{2}}=\sum_{i=1,2} \sigma_i\lambda_{i}\left(1+\lambda_{i}^{2}\right)^{\frac{1}{4}}E\left(\frac{\pi}{2},k_{i}\right)+\sum_{i=1,2} \sigma_i \frac{(1+\lambda_i^2)^{1/4} }{ 2 ( \lambda_{i}+\sqrt{1+\lambda_{i}^{2} })} F(\frac{\pi}{2}, k_i ) $$
and
$$ \sqrt{\Delta}\sum_{i=1,2}\biggl[-\left(1+\lambda_{i}^{2}\right)^{\frac{1}{4}}E\left(\frac{\pi}{2},k_{i}\right)+\frac{F\left(\frac{\pi}{2},k_{i}\right)}{2\left(1+\lambda_{i}^{2}\right)^{\frac{1}{4}}}\left( \lambda_{i} + \frac{1}{\sqrt{1+\lambda_{i}^{2}}+\lambda_{i}} \right)\biggr]=\sum_{i=1,2}\biggl[-\left(1+\widetilde{\lambda}_{i}^{2}\right)^{\frac{1}{4}}E\left(\frac{\pi}{2},\widetilde{k}_{i}\right)+\frac{F\left(\frac{\pi}{2},\widetilde{k}_{i}\right)}{2\left(1+\widetilde{\lambda}_{i}^{2}\right)^{\frac{1}{4}}}\left(\widetilde{\lambda}_{i}+\frac{1}{\sqrt{1+\widetilde{\lambda}_{i}^{2}}+\widetilde{\lambda}_{i}} \right)\biggr] $$
where $\sigma_1=1$, $\sigma_2=-1$,
$k_{i}^{2}=\frac{\sqrt{1+\lambda_{i}^{2}}+\lambda_{i}}{2\sqrt{\left(1+\lambda_{i}^{2}\right)}}$ for $i=1,2$ with $\lambda_{1}\equiv \frac{\mu-\Delta_{0}}{\Delta}$, $\lambda_{2}\equiv \frac{-\Delta_{0}-\mu}{\Delta} $, $\widetilde{k}_{i}^{2}=\frac{\sqrt{1+\widetilde{\lambda}_{i}^{2}}+\widetilde{\lambda}_{i}}{2\sqrt{\left(1+\widetilde{\lambda}_{i}^{2}\right)}}$ for $i=1,2$ with $\widetilde{\lambda}_{1}\equiv \mu$, $\widetilde{\lambda}_{2}\equiv -\mu $. $F(\frac{\pi}{2},k)$ and $E(\frac{\pi}{2},k)$ are the complete elliptic integral of the first and the second kind respectively which in Mathematica are EllipticK[] and EllipticE[]. In the code below $\mu$, $\Delta_0$, $\Delta$, $\varepsilon$ are represented by u, o, d, e.
By first converting those two equations in the form of differential equations and then using NDSolveValue to obtain plots of $\Delta$ and $\mu$ as a function of $\Delta_0$ at $\varepsilon=0.001$, I encounter the error:
lambda1[u_, o_, d_] := (u - o)/d;
lambda2[u_, o_, d_] := (-u - o)/d;
k1sq[u_, o_, d_] := (Sqrt[1 + ((lambda1[u, o, d])^2)] + lambda1[u, o, d])/(2*(Sqrt[1 + ((lambda1[u, o, d])^2)]));
k2sq[u_, o_, d_] := (Sqrt[1 + ((lambda2[u, o, d])^2)] + lambda2[u, o, d])/(2*(Sqrt[1 + ((lambda2[u, o, d])^2)]));
a1[u_, o_, d_] := (1 + (lambda1[u, o, d])^2)^(1/4);
a2[u_, o_, d_] := (1 + (lambda2[u, o, d])^2)^(1/4);
b1[u_, o_, d_] := Sqrt[1 + (lambda1[u, o, d])^2] + lambda1[u, o, d];
b2[u_, o_, d_] := Sqrt[1 + (lambda2[u, o, d])^2] + lambda2[u, o, d];
equat1[u_, o_, d_, e_] =
(2/3)*(((e)/d)^(3/2)) == ((lambda1[u, o, d]*a1[u, o, d]*EllipticE[k1sq[u, o, d]])
-((lambda1[u, o, d]*EllipticK[k1sq[u, o, d]])/(2*a1[u, o, d]*b1[u, o, d]))+
((EllipticK[k1sq[u, o, d]])/(2*a1[u, o, d])) -
(lambda2[u, o, d]*a2[u, o, d]* EllipticE[k2sq[u, o, d]])
+ ((lambda2[u, o, d]*EllipticK[k2sq[u, o, d]])/(2*a2[u, o, d]*b2[u, o, d]))-
((EllipticK[k2sq[u, o, d]])/(2*a2[u, o, d])));
equat2[u_, o_, d_] =
Sqrt[d]*(((lambda1[u, o, d])*
EllipticK[k1sq[u, o, d]]/(2*a1[u, o, d])) - (a1[u, o, d]*
EllipticE[
k1sq[u, o, d]]) + ((EllipticK[k1sq[u, o, d]])/(2*a1[u, o, d]*
b1[u, o, d])) + ((lambda2[u, o, d])*
EllipticK[k2sq[u, o, d]]/(2*a2[u, o, d])) - (a2[u, o, d]*
EllipticE[
k2sq[u, o, d]]) + ((EllipticK[k2sq[u, o, d]])/(2*a2[u, o, d]*
b2[u, o, d]))) == ((lambda1[u, 0, 1])*
EllipticK[k1sq[u, 0, 1]]/(2*a1[u, 0, 1])) - (a1[u, 0, 1]*
EllipticE[
k1sq[u, 0, 1]]) + ((EllipticK[k1sq[u, 0, 1]])/(2*a1[u, 0, 1]*
b1[u, 0, 1])) + ((lambda2[u, 0, 1])*
EllipticK[k2sq[u, 0, 1]]/(2*a2[u, 0, 1])) - (a2[u, 0, 1]*
EllipticE[
k2sq[u, 0, 1]]) + ((EllipticK[k2sq[u, 0, 1]])/(2*a2[u, 0, 1]*
b2[u, 0, 1]));
{dk, uk} = {d, u} /. FindRoot[{equat2[u, o, d], equat1[u, o, d, .001]} /. o -> 0, {{d, .5}, {u, .1}}];
{dsolk, usolk} =
NDSolveValue[{D[equat2[u, o, d] /. {d -> d[o], u -> u[o]}, o],
D[equat1[u, o, d, .001] /. {d -> d[o], u -> u[o]}, o],
u[0] == uk, d[0] == dk}, {d, u}, {o, 0, 1},
Method -> {"EquationSimplification" -> "Residual"}];
Plot[{dsolk[t], usolk[t]}, {t, 0, 1}, PlotRange -> {-.1, 1.1}, Frame -> True]
NDSolveValue::mxst: Maximum number of 100000 steps reached at the point o == 0.6097033678408679`.
I expect the blue curve ($\Delta$) to continue levelling off and the yellow curve ($\mu$) to continue increasing. Increasing MaxSteps->Infinity probably helps but the run-time took too long to evaluate and I gave up. Eventually, I will be obtaining these plots for different values of $\varepsilon$ and so increasing the run-time by increasing MaxSteps does not seem to be the way to go.
Using FindRoot directly gives somwhat the plot I want and the run-time is relatively shorter, but fails around the sharp turn at $\Delta_0=0.29$ and throws the error:
Table[Values@
FindRoot[{equat2[u, o, d],
equat1[u, o, d, .001]}, {{d, .5}, {u, .1}}], {o, 0, 1, .001}];
ListPlot[Transpose@%, DataRange -> {0, 1}, Joined -> True,
AxesLabel -> {o, "d, u"},
LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]
FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.
