Definitions

SetDirectory[NotebookDirectory[]];

SetAttributes[untilSuccess,HoldAll];
untilSuccess[e_]:=Quiet[Check[e,untilSuccess[e]]]

μm[f_]:=(Sqrt[4f''[1]])

Eopt[f_,μ_]:=((μ(f[1] f'[1]-((f'[1])^2)-f[1] f''[1]))/(2 f'[1] f''[1]))

Eth[f_]:=Eopt[f,μm[f]]

Eμopt[f_,μ_]:=-((μ f'[1])/(f'[1]+f''[1]))

𝒟[f_][μ_]:=With[{μm=(Sqrt[4f''[1]])},(1/2) +(((μ )^2)/((μm)^2))+Log[(1/2)μm]-UnitStep[(μ^2)-((μm)^2)](Log[RealAbs[(μ/(μm))]-(Sqrt[((μ^2)/((μm)^2))-1])]+ RealAbs[(μ/(μm))](Sqrt[((μ^2)/((μm)^2))-1]))]
d𝒟[f_][μ_]:=With[{μm=(Sqrt[4f''[1]])},(2/(μm))((μ/(μm))-UnitStep[(μ^2)-((μm)^2)]RealSign[μ] (Sqrt[((μ^2)/((μm)^2))-1]))]

Σa[f_][E_,μ_]:=𝒟[f][μ]-(1/2)-(1/2)Log[f'[1]]-(1/(2f[1]))((E^2)+((((f[1] μ+f'[1] E))^2)/(f[1] (f'[1]+f''[1])-((f'[1])^2))))

ΣaSaddles[f_][EE_,μ_]=Simplify[Σa[f][EE,μ],Assumptions->{(μ^2)<((μm[f])^2),μ>0}];
ΣaMinima[f_][EE_,μ_]=Simplify[Σa[f][EE,μ],Assumptions->{(μ^2)>((μm[f])^2),μ>0,f''[1]>0}];

𝒮tmp1[f_][E_,μ_][q1_,x_][β_,rd_,r1_,dd_,d1_]:=With[{Δx=1-x},
β E-rd μ+(1/2)((β^2)(f[1]-Δx f[q1])+(2 β rd-dd)f'[1]-Δx(2 β r1-d1)f'[q1]+(rd^2)f''[1]-Δx (r1^2)f''[q1]+(1/x)Log[(((rd-Δx r1))^2)+dd(1-Δx q1)-Δx d1(1-Δx q1)]-((Δx)/x)Log[(((rd-r1))^2)+(dd-d1)(1-q1)])
]

𝒮tmp2[f_][EE_,μ_][q1_,x_][β_,rd_,r1_]=Simplify[𝒮tmp1[f][EE,μ][q1,x][β,rd,r1,dd,d1]/.First[Solve[0==D[𝒮tmp1[f][EE,μ][q1,x][β,rd,r1,dd,d1],{{dd,d1}}],{dd,d1}]],Assumptions->{((1-q1 (1-x))/(f'[1]-(1-x) f'[q1]))>0,((1-q1)/(f'[1]-f'[q1]))>0,0<x<1,0<q1<1}];

𝒮[f_][EE_,μ_][q1_,x_]=𝒟[f][μ]+Simplify[𝒮tmp2[f][EE,μ][q1,x][β,rd,r1]/.First[Solve[0==D[𝒮tmp2[f][EE,μ][q1,x][β,rd,r1],{{β,rd,r1}}],{β,rd,r1}]]];

u[f_]:=f[1](f'[1]+f''[1])-((f'[1])^2)
v[f_]:=f'[1](f''[1]+f'''[1])-((f''[1])^2)
w[f_]:=2((f''[1])^2)+f'[1](f'''[1]-2f''[1])
y[f_]:=f'[1](f'[1]-f[1])+f''[1]f[1]
z[f_]:=f[1](f''[1]-f'[1])+((f'[1])^2)

μ0[f_][E_]:=((-1)/(z[f]))(2E f'[1]f''[1]+(Sqrt[2f''[1] u[f](Log[((f''[1])/(f'[1]))]z[f]-(E^2)(f''[1]-f'[1]))]))

Emax[f_]=EE/.Solve[μ0[f][EE]==μm[f],EE][[1]]//Simplify

q1Exp={f[q1]->f[1]+f'[1](q1-1)+(1/2)f''[1](((q1-1))^2),Derivative[n_][f][q1]:>Derivative[n][f][1]+Derivative[n+1][f][1](q1-1)+(1/2)Derivative[n+2][f][1](((q1-1))^2)};

d𝒮dq1dx=Simplify[Limit[Simplify[D[D[𝒮[f][ϵ,μ][q1,x],q1],x]/.x->1]/.q1Exp,q1->1,Direction->"FromBelow"]];

d𝒮Coeff=FullSimplify[CoefficientList[d𝒮dq1dx/.μ->((-1)/(z[f]))(2ϵ f'[1]f''[1]+y),{ϵ,y}]]

ΣaCoeff=FullSimplify[CoefficientList[(2 f[1] f'[1]-2 ((f'[1])^2)-2 f[1] f''[1])ΣaSaddles[f][ϵ,((-1)/(z[f]))(2ϵ f'[1]f''[1]+y)],{ϵ,y}]]

letterRules={d->d𝒮Coeff[[1,1]],c->d𝒮Coeff[[2,2]],b->d𝒮Coeff[[3,1]],a->d𝒮Coeff[[1,3]],g->-ΣaCoeff[[1,1]],e->ΣaCoeff[[3,1]],h->ΣaCoeff[[1,3]]};

solLetters=Solve[{0==a (y^2)+b (ϵ^2)+c ϵ y+d,0==h (y^2)+e (ϵ^2)-g},{y,ϵ}]//Simplify;

solLetters

letterRules

E1RSBp[f_]:=Module[{
a=(-8 f[1] (f[1]-f'[1]) (f'[1]-f''[1]) ((f''[1])^2)+(3 (((f[1]-f'[1]))^2) ((f'[1])^2)+2 f[1] f'[1] (-f[1]+f'[1]) f''[1]+3 ((f[1])^2) ((f''[1])^2)) Derivative[3][f][1])/(16 ((((((f[1]-f'[1]))^2) ((f'[1])^2) f''[1]-((f[1])^2) ((f''[1])^3)))^2)),
b=((f'[1] (2 ((f''[1])^2)+f'[1] (-2 f''[1]+Derivative[3][f][1])))/(4 (((((f'[1])^2)+f[1] (-f'[1]+f''[1])))^2))),
c=((2 ((f''[1])^2)+f'[1] (-2 f''[1]+Derivative[3][f][1]))/(4 f''[1] (((((f'[1])^2)+f[1] (-f'[1]+f''[1])))^2))),
d=(1/4) (2-((2 f''[1])/(f'[1]))-((Derivative[3][f][1])/(f''[1]))),
e=-f'[1]+f''[1],
g=-((Log[f'[1]]-Log[f''[1]]) (((f'[1])^2)+f[1] (-f'[1]+f''[1]))),
h=(1/(2 f''[1] (-((f'[1])^2)+f[1] (f'[1]+f''[1]))))
},
{-(1/(2 (Sqrt[2]) c e (d e+b g)))(Sqrt[(((c^2) e g-2 a e (d e+b g)+2 b d e h+2 (b^2) g h-(Sqrt[(c^2) (e^2) ((c^2) (g^2)-4 a g (d e+b g)-4 d (d e+b g) h)]))/((a^2) (e^2)-2 a b e h+h ((c^2) e+(b^2) h)))]) ((c^2) e g+(Sqrt[(c^2) (e^2) ((c^2) (g^2)-4 a g (d e+b g)-4 d (d e+b g) h)])),
(1/(2 (Sqrt[2]) c e (d e+b g)))(-(c^2) e g+(Sqrt[(c^2) (e^2) ((c^2) (g^2)-4 a g (d e+b g)-4 d (d e+b g) h)])) (Sqrt[(((c^2) e g-2 a e (d e+b g)+2 b d e h+2 (b^2) g h+(Sqrt[(c^2) (e^2) ((c^2) (g^2)-4 a g (d e+b g)-4 d (d e+b g) h)]))/((a^2) (e^2)-2 a b e h+h ((c^2) e+(b^2) h)))])}
]

G[f_]:=-2(f''[1]-f'[1])u[f]w[f]-2((Log[((f''[1])/(f'[1]))])^2)((f'[1])^2)f''[1]v[f]-f'[1]Log[((f''[1])/(f'[1]))](4(f'[1]-2f[1])(f''[1]-f'[1])((f''[1])^2)-(-3(f'[1]-f[1])((f'[1])^2)+f'[1](f'[1]-2f[1])f''[1]+3f[1]((f''[1])^2))f'''[1])

h11Auff[f_,x_]:=((x f'[1]-f'[x])/(f'[1]-f'[x]))+(((((f'[1]x-f'[x]))^2)f[x])/(f'[1]f'[x](f'[1]-f'[x])x(1-x)))-Log[((f'[1]-f'[x])/(f'[1](1-x)))]
h12Auff[f_,x_]:=(1-x)(f'[1]-f'[x])(((f'[1]f'[x]x)/(((f'[1]x-f'[x]))^2))Log[((f'[1]x)/(f'[x]))]-((f'[x])/(f'[1]x-f'[x])))+f'[x](1-x)-1+f[x]
Ψ[λ_,p_,s_]:=With[{f=Function[q,(λ (q^p)+(1-λ)(q^s))]},-((f'[1]-1)/(f''[1]))-Log[((f''[1])/(f'[1]))]+1-(2/(f'[1]))+((f''[1])/((f'[1])^2))]

cleanContourPlot[cp_Graphics]:=Module[{points,groups,regions,lines},groups=Cases[cp,{style__,g_GraphicsGroup}:>{{style},g},Infinity];
points=First@Cases[cp,GraphicsComplex[pts_,___]:>pts,Infinity];
regions=Table[Module[{group,style,polys,edges,cover,graph},{style,group}=g;
polys=Join@@Cases[group,Polygon[pt_,___]:>pt,Infinity];
edges=Join@@(Partition[#,2,1,1]&/@polys);
cover=Cases[Tally[Sort/@edges],{e_,1}:>e];
graph=Graph[UndirectedEdge@@@cover];
{Sequence@@style,FilledCurve[List/@Line/@First/@Map[First,FindEulerianCycle/@(Subgraph[graph,#]&)/@ConnectedComponents[graph],{3}]]}],{g,groups}];
lines=Cases[cp,_Tooltip,Infinity];
Graphics[GraphicsComplex[points,{regions,lines}],Sequence@@Options[cp]]]

cleanContourPlot[Legended[cp_Graphics,rest___]]:=Legended[cleanContourPlot[cp],rest]
cleanRegionPlot[cp_Graphics]:=Module[{points,groups,regions,lines},groups=Cases[cp,{style__,g_GraphicsGroup}:>{{style},g},Infinity];
points=First@Cases[cp,GraphicsComplex[pts_,___]:>pts,Infinity];
regions=Table[Module[{group,style,polys,edges,cover,graph},{style,group}=g;
polys=Join@@Cases[group,Polygon[pt_,___]:>pt,Infinity];
edges=Join@@(Partition[#,2,1,1]&/@polys);
cover=Cases[Tally[Sort/@edges],{e_,1}:>e];
graph=Graph[UndirectedEdge@@@cover];
{Sequence@@style,FilledCurve[List/@Line/@First/@Map[First,FindEulerianCycle/@(Subgraph[graph,#]&)/@ConnectedComponents[graph],{3}]]}],{g,groups}];
lines=Cases[cp,{__,_Line},Infinity];(*only change*)Graphics[GraphicsComplex[points,{regions,lines}],Sequence@@Options[cp]]];

eqZero[f_,k_]:=-(1/2)Rd (OverHat[β])f'[1]-(1/2)(((OverHat[β])^2)(f[1]+Sum[(x[i]-x[i+1])f[q[i]],{i,0,k}])+(1/(x[1]))Log[1+(OverHat[β])(Rd^(-1))(1+Sum[(x[i]-x[i+1]) q[i],{i,0,k}])]+Sum[((1/(x[j+1]))-(1/(x[j])))Log[1+(OverHat[β])(Rd^(-1))(1+Sum[(x[i]-x[i+1]) q[i],{i,j+1,k}]-x[j+1]q[j])],{j,1,k}])/.{x[0]->0,q[0]->0,x[k+1]->1}

nMaximize[f_,k_,n_]:=
First@Reverse@SortBy[
ParallelTable[untilSuccess@FindMaximum[Prepend[Flatten@Join[{Rd>0,(OverHat[β])>0,0<q[1],0<x[1],q[k]<1,x[k]<1},Table[{0<q[i+1]-q[i],0<x[i+1]-x[i]},{i,1,k-1}]],(1/(OverHat[β]))eqZero[f,k]],Thread[{Join[x/@Range[k],q/@Range[k],{(OverHat[β]),Rd}],RandomReal[{0,1},2k+2,WorkingPrecision->50]}],WorkingPrecision->50],{n}],
First]

fps[p_,s_][λ_]:=Function[q,(1/2)(λ (q^p)+(1-λ)(q^s))]

Figures

Phases figure

λAuffStart=FindRoot[Evaluate[{Ψ[λ,3,s],D[Ψ[λ,3,s],λ]}],{{s,12.2},{λ,0.86}},WorkingPrecision->20]

λ12=Prepend[Rest[FoldList[Function[{sol,ss},Prepend[FindRoot[{h11Auff[Function[q,(λ (q^3)+(1-λ)(q^ss))],x],h12Auff[Function[q,(λ (q^3)+(1-λ)(q^ss))],x]},Rest[sol]/.Rule->List,MaxIterations->1000,WorkingPrecision->20],s->ss]],{s->(s/.λAuffStart),λ->(0.8901),x->0.9664},s+Rest@Range[0,8,1/100]/.λAuffStart]],λAuffStart];

λ21=Prepend[Rest[FoldList[Function[{sol,ss},Prepend[FindRoot[Ψ[λ,3,ss],Rest[sol]/.Rule->List,MaxIterations->1000,WorkingPrecision->10],s->ss]],{s->(s/.λAuffStart),λ->(λ+0.01/.λAuffStart)},s+Rest@Range[0,8,1/100]/.λAuffStart]],λAuffStart];

rsbComplexityStart=FindRoot[D[G[Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))]],λ]/.λ->0,{s,4.5},WorkingPrecision->20]

rsbComplexityBoundary=Prepend[Rest[FoldList[Function[{sol,ss},Prepend[FindRoot[G[Function[q,(1/2)(λ (q^3)+(1-λ)(q^ss))]],Rest[sol]/.Rule->List,WorkingPrecision->15],s->ss]],Append[rsbComplexityStart,λ->(10^(-1))],s+Rest[Range[0,16,1/100]]/.rsbComplexityStart]],Append[rsbComplexityStart,λ->0]];

plot3s=Show[
ListPlot[{{s,λ}/.rsbComplexityBoundary,{s,λ}/.λ12,{s,λ}/.λ21,{{8,0},{8,1.1},{20,1.1}}},Joined->True,Filling->{1->0,2->{3},4->0},PlotStyle->{ColorData[97][1],ColorData[97][3],ColorData[97][3],ColorData[97][2]},AspectRatio->1/GoldenRatio,PlotRange->{{3,17},{0,1}},FrameLabel->{s,λ},RotateLabel->False,LabelStyle->Directive[Black,FontFamily->"Times",FontSize->10,Filling->{{1}->0}],FrameStyle->Black,Frame->True],
Epilog->{Text[Style["RSB\n complexity",Directive[ColorData[97][1],FontFamily->"Times",FontSize->10]],{6.2,0.3}],
Text[Style["Non\[Hyphen]convex f",Directive[ColorData[97][2],FontFamily->"Times",FontSize->10]],{12.5,0.4}],
Text[Style["RSB T = 0",Directive[ColorData[97][3],FontFamily->"Times",FontSize->10]],{15,0.85}]},ImageSize->270,RotateLabel->False
]

Export["figs/phases_34.pdf",plot3s];

Complexity figure

f35=Function[q,(1/2)((1/2) (q^3)+(1-(1/2))(q^5))];

Eini=E1RSBp[f35][]-(10^(-3));

zeroComplexityEqs=Append[D[𝒮[f35][ϵ,μ][q1,x],{{q1,x}}],𝒮[f35][ϵ,μ][q1,x]];

μ0[f35][Eini]//N

iniComplexityZero=Select[Table[FindRoot[zeroComplexityEqs/.ϵ->Eini,{{μ,μ0[f35][Eini]},{q1,1-0.01RandomReal[],0,1},{x,1-0.01RandomReal[],0,1}},WorkingPrecision->15],{100}],Norm[zeroComplexityEqs/.ϵ->Eini]<(10^(-4))&&x<1/.#&];

Length[iniComplexityZero]

iniComplexityZero2=FindRoot[zeroComplexityEqs/.ϵ->Eini,(#1->Re[#2]&@@@RandomChoice[iniComplexityZero])/.Rule->List,WorkingPrecision->30]

solsComplexityZero=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->f35,ϵ->ϵϵ},FindRoot[zeroComplexityEqs/.ϵ->ϵϵ,Drop[sol,2]/.Rule->List,WorkingPrecision->20]]],
Join[{f->f35,ϵ->Eini},iniComplexityZero2],Eini-Rest@Range[0,1,(10^(-4))],0<=x<=1&&0<=q1<=1&&Norm[zeroComplexityEqs]<(10^(-6))//.#&]];

transitionEqs=Simplify[{D[𝒮[f35][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[f35][ϵ,μ][q1,x],q1],{x,1,1}]}];

transitionIni=FindRoot[transitionEqs/.ϵ->Eini,{{μ,μ0[f35][Eini]},{q1,1-0.1RandomReal[]}},WorkingPrecision->20]

solsTransition=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->f35,ϵ->ϵϵ,x->1},FindRoot[transitionEqs/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->f35,ϵ->Eini,x->1},transitionIni],Eini-Range[0,2,(10^(-4))],0<=q1<=1&&𝒮[f35][ϵ,μ][q1,x]>=0/.#&]];

transitionEnd=FindRoot[Append[transitionEqs,𝒮[f35][ϵ,μ][q1,1]],Delete[solsTransition[[-1]],{{1},{3}}]/.Rule->List,WorkingPrecision->50,MaxIterations->500]

(ϵ/.transitionEnd)-E1RSBp[f35][-1]

zoomPlot=Show[
ListPlot[{{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.solsTransition,{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.solsComplexityZero,{{-1.3,0},{0.1,0}}},Joined->True,PlotRange->{{-1.25,-0.85},{-0.25,1.15}},Epilog->{Red,Point[{#,0}]&/@{E1RSBp[f35][],E1RSBp[f35][-1]}},PlotStyle->{Red,Black,{Black,Dashed}},Filling->{2->{1}},Prolog->{Red,Opacity[0.3],Polygon[Join[{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.solsTransition,{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.Reverse@solsComplexityZero]]},Frame->True,FrameStyle->Black,LabelStyle->Directive[Black,FontFamily->"Times",FontSize->10],FrameLabel->{None,"(μ - (Subscript[μ, 0])) × (10^3)"},ImageSize->140,FrameTicks->{{{0,1},None},{{-1.2,-1.1,-1,-0.9},None}},Axes->None]
]

nonZoomPlot=ListPlot[{{ϵ,μ}/.solsTransition,{ϵ,μ}/.solsComplexityZero},Joined->True,PlotStyle->Red,Epilog->{Red,Point[{#,μ0[f35][#]}]&/@{E1RSBp[f35][],E1RSBp[f35][-1]}}]

rp=RegionPlot[Evaluate[Σa[f35][ϵ,μ]>0&&μ>μm[f35]],{ϵ,-1.5,-1.1},{μ,4.5,5.5},PlotPoints->100,PlotStyle->{Opacity[0.5],Gray},BoundaryStyle->None]

plot35=Show[ContourPlot[0==Σa[f35][ϵ,μ],{ϵ,-1.3,-0.8},{μ,2,5.5},ContourStyle->{{Black,Dashed}}],rp,nonZoomPlot,ImageSize->270,Prolog->Inset[zoomPlot,{-1.13,2.7}],LabelStyle->Directive[Black,FontFamily->"Times",FontSize->10],FrameStyle->Black,FrameLabel->{"E",μ},RotateLabel->False,Epilog->{Red,Point[{#,μ0[f35][#]}]&/@{E1RSBp[f35][],E1RSBp[f35][-1]}}]

Export["figs/complexity_35.pdf",plot35];

Ranges figure

fTest=Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]/.λ->650/1000;

Eini=E1RSBp[fTest][[1]]+1/10000;

transitionEqs=Simplify[{D[𝒮[fTest][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[fTest][ϵ,μ][q1,x],q1],{x,1,1}]}];

transitionIni=FindRoot[transitionEqs/.ϵ->Eini,{{μ,μ0[fTest][Eini]},{q1,1-0.1RandomReal[]}},WorkingPrecision->20]

solsTransition2=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest,ϵ->ϵϵ,x->1},FindRoot[transitionEqs/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest,ϵ->Eini,x->1},transitionIni],Eini+Range[0,2,(10^(-3))/2],0<=q1<=1&&𝒮[fTest][ϵ,μ][q1,x]>=0/.#&]];

solsTransition=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest,ϵ->ϵϵ,x->1},FindRoot[transitionEqs/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest,ϵ->Eini,x->1},transitionIni],Eini-Range[0,2,(10^(-4))],0<=q1<=1&&𝒮[fTest][ϵ,μ][q1,x]>=0/.#&]];

Show[
ListPlot[{{ϵ,μ}/.solsTransition2,{ϵ,μ}/.solsTransition},Joined->True,PlotRange->All],
Plot[μ0[fTest][ϵ],{ϵ,-2,2}]
]

transitionEnd=FindRoot[Append[transitionEqs,𝒮[fTest][ϵ,μ][q1,1]],Delete[solsTransition[[-1]],{{1},{3}}]/.Rule->List,WorkingPrecision->50,MaxIterations->500]

E1RSBp[fTest][[1]]-ϵ/.transitionEnd

fλ=Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))];

transitionEqsλ=Simplify[{D[𝒮[fλ][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[fλ][ϵ,μ][q1,x],q1],{x,1,1}]}];

solsEnd1=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],Prepend[transitionEnd,λ->65/100],Rest[Range[65/100,9/10,(10^(-2))/5]]];

solsEnd2=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],solsEnd1[[-1]],Rest[Range[9/10,945/1000,(10^(-3))/5]]];

solsEnd3=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],solsEnd2[[-1]],Rest[Range[945/1000,960/1000,(10^(-5))]]];

solsEnd4=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],Prepend[transitionEnd,λ->65/100],Rest[Range[65/100,25/100,-(10^(-2))/5]]];

solsEnd5=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],solsEnd4[[-1]],Rest[Range[25/100,24/100,-(10^(-4))]]];

solsEndAll=Join[Reverse[solsEnd5],Reverse[solsEnd4],solsEnd1,solsEnd2,solsEnd3];

ListPlot[{λ,μ-μm[fλ]}/.solsEndAll,PlotRange->All,Joined->True]

(s|->Show[
Plot[Evaluate[Append[E1RSBp[f],Emax[f]]/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))]],{λ,0,1},PlotStyle->{Red,Red,Black},Filling->{1->{2}}]
])/@{5,8,11,14}

rangePlotsLog=(s|->Show[
LogPlot[Evaluate[Max[#,(10^(-10))]&/@(E1RSBp[f]-Emax[f]/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))])],{λ,0,1},PlotStyle->{Red,Red},Filling->{1->{2}},PlotRange->{(10^(-6)),1.5}]
])/@{5,8,11,14}

rangePlotsLogμ=(s|->Show[
LogPlot[Evaluate[Max[#,(10^(-10))]&/@({μm[f]-μ0[f][E1RSBp[f][[2]]],(μm[f]-μ0[f][E1RSBp[f][[1]]])If[E1RSBp[f][[1]]<Emax[f],0,1]}/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))])],{λ,0,1},PlotStyle->{Red,Red},Filling->{1->{2}},PlotRange->{(10^(-6)),20}]
])/@{5,8,11,14}

Show[rangePlotsLogμ[[4]],ListLogPlot[{λ,Max[μm[fλ]-μ,(10^(-10))]}/.solsEndAll,PlotRange->All,Joined->True]]

Show[rangePlotsLog[[4]],ListLogPlot[{λ,Max[ϵ-Emax[fλ],(10^(-10))]}/.solsEndAll,PlotRange->All,Joined->True]]

fTest=Function[q,(1/2)(λ (q^3)+(1-λ)(q^7))]/.λ->200/1000;

Eini=E1RSBp[fTest][]+1/100000;

transitionEqs=Simplify[{D[𝒮[fTest][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[fTest][ϵ,μ][q1,x],q1],{x,1,1}]}];

E1RSBp[fTest][]//N

q1Exp={f[q1]->f[1]+f'[1](q1-1)+(1/2)f''[1](((q1-1))^2),Derivative[n_][f][q1]:>Derivative[n][f][1]+Derivative[n+1][f][1](q1-1)+(1/2)Derivative[n+2][f][1](((q1-1))^2)}

μ/.solLetters[[1]]//FullSimplify

ϵ/.solLetters[[1]]//FullSimplify

ϵ/.solLetters/.letterRules/.f->fTest//N

E1RSBp[fTest][]//N

μ0[fTest][E1RSBp[fTest][]]//N

μ/.solLetters/.letterRules/.f->fTest//N

testest=Solve[0==ΣaSaddles[f][ϵ,μ],μ]//FullSimplify[#,Assumptions->{-((f'[1])^2)+f[1] (f'[1]+f''[1])>0,f'[1]>0,f''[1]>0}]&

fTest=Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]/.λ->1/20

Eini=E1RSBp[fTest][[2]]-1/1000;

transitionIni=FindRoot[transitionEqs/.ϵ->Eini,{{μ,μ0[fTest][Eini]},{q1,1-0.01RandomReal[]}},WorkingPrecision->20]

solsTransition=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest,ϵ->ϵϵ,x->1},FindRoot[transitionEqs/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest,ϵ->Eini,x->1},transitionIni],Eini-Range[0,2,(10^(-4))],0<=q1<=1&&𝒮[fTest][ϵ,μ][q1,x]>=0/.#&]];

transitionEnd=FindRoot[Append[transitionEqs,𝒮[fTest][ϵ,μ][q1,1]],Delete[solsTransition[[-1]],{{1},{3}}]/.Rule->List,WorkingPrecision->50,MaxIterations->500]

ϵ-Emax[fTest]/.transitionEnd

Eini=E1RSBp[fTest][-1]-1/1000;

transitionOtherIni=FindRoot[transitionEqs/.ϵ->Eini,{{μ,μ0[fTest][Eini]},{q1,1-0.1RandomReal[]}},WorkingPrecision->20]

solsTransition=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest,ϵ->ϵϵ,x->1},FindRoot[transitionEqs/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest,ϵ->Eini,x->1},transitionIni],Eini-Range[0,2,(10^(-4))],0<=q1<=1&&𝒮[fTest][ϵ,μ][q1,x]>=0/.#&]];