Integrate out and define the vacuum for the fundamental --- (Subscript[m, q|UV])=0.

Defining Constants

rmax=10;(*UV cuttoff*)
nc=5;(*Number of colours*) (*working in SU(2,3,4,5)*)
c2g=nc;(*for the β-function*) (*for this theory, SU(N)*)
nfr1=16.9;(*fundamental flavours*) (*varying 1 to N_f(max)=6.5,10,13.5,17 for SU(2,3,4,5)*)
dr1=nc;(*dimension of the fundamental*)
tr1=(1/2);(*half the Dynkin index of the fundamental*)
c2r1=(((nc^2)-1)/(2nc));(*quadatic Casimir of the fundamental*)
nfr2=(1/2);(*adjoint flavours*)
dr2=(nc^2)-1;(*dimension of the adjoint*)
tr2=nc;(*half the Dynkin index of the adjoint*)
c2r2=nc;(*quadatic Casimir of the adjoint*)

Symmetric

Nfp=nfr2;(*number of flavours used to calculate fπ*)
g5sq:=((12 Pi^2)/(dr2 nfr2));
α1:=(π/(3.*c2r2)); (*Value of α[1] compare to output at end to tune for δ*)
B0:=(1/(6 π))(11 c2g -4 tr1 nfr1 -4 tr2 nfr2) ;(*one-loop*) (*doubled factors to 4 because Dirac fermions*)
B1:=(1/(24 (π^2)))(34(c2g^2)- (20c2g + 12 c2r1)tr1 nfr1- (20c2g + 12c2r2)tr2 nfr2 );(*two-loop*) (*we have doubled the factors to 20 and 12*)
eqalp1:=alp1'[mu1] + (B0 alp1[mu1]^2  + B1 alp1[mu1]^3)==0
sol1 := NDSolve [{ eqalp1, alp1[0]==α1}, alp1, {mu1,-100, (10^24)},Method->{"StiffnessSwitching"}];
alph1[mu1_]:= First[Evaluate[alp1[mu1]/. sol1]]
αv1[mu_]:=alph1[mu]
γ1[mu1_]:=(3/(2π)) *c2r2*αv1[mu1] 
δ1[mu1_]:=-2γ1[mu1]
δ1[0]
del1[mu1_]=-2γ1[mu1];

mρ=Vmass1;
plotdel1v1=Plot[del1[μ+Log[mρ]],{μ,-100,100},PlotStyle->{Thickness[0.009],Blue},Frame->True,FrameLabel->{Style["ln (μ/(Subscript[m, ρ, 15]))",Black,17],Style["((Δm)^2)",Black,17]},FrameTicksStyle->Directive[Black,17],PlotTheme->"Scientific",FrameStyle->{{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black}},BaseStyle->{17,FontFamily->"Times New Roman"},AxesStyle->{Thick,Black},PerformanceGoal->"Quality",PlotRange->{{-60,60},{-4,0}}]

mρ=Vmass1;
muv=0.0;(*UV mass (IF using IR mass comment this and the line to calculate mir later)*)
guess1=3.*(10^(-7));
ode:=D[(r^3)L'[r],r]-δ1[Log[Sqrt[(r^2)+((L[r])^2)] ]]r L[r]==0
solm[v_]:= NDSolve[{ode,L[v ]==v,
L'[v]==0},L,{r,v,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First
LinfBm[ v_?NumericQ ]:=(L[rmax]/.solm[v])-muv
(*Find mir from muv (comment out if using mir*)
mir1 = m/.FindRoot[LinfBm[m],{m,guess1}]
lf1[ρ_]:=Evaluate[L[ρ]/.solm[ mir1]]
plot1=Show[Plot[((lf1[r*mρ])/(mρ)),{r,(mir1/(mρ)),rmax},PlotStyle->{Thickness[0.009],Blue},Frame->True,FrameLabel->{Style["(ρ/(Subscript[m, ρ, 15]))",Black,Bold,17],Style["((L ((Subscript[ρm, ρ, 15])))/(Subscript[m, ρ, 15]))",Black,Bold,17]},FrameTicksStyle->Directive[Black,17],PlotTheme->"Scientific",FrameStyle->{{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black}},BaseStyle->{17,FontFamily->"Times New Roman"},AxesStyle->{Thick,Black},PerformanceGoal->"Quality", PlotRange->{{0,rmax},{-0.01,0.4}},PlotLegends->"Blue: Vacuum for the symmetric"],Plot[x,{x,0,6}, PlotStyle->{Thick,Black}]]
(*The above is the "embedding"*)
(*To get the point that we set the wall*)
(*It should be the IR-mass, but check for better or for worse*)
{r,lf1[r]}/.FindRoot[lf1[r]==r,{r,guess1}]
(*Now, we can write L(ρ=(Subscript[ρ, IR]))=(Subscript[ρ, IR]) and L'(ρ=(Subscript[ρ, IR]))=0*)
(*Create a table of elements and test the p|ositivity of the plot*)
table1=Table[lf1[r],{r,mir1,rmax,(10^(-2))}];
AllTrue[table1,#>0&]
lf1[r]/.r:>10//SetPrecision[#,2]&
L1[ρ_]=lf1[ρ];

The  equations  of  motion - Symmetric

Vode:=D[(r^3)V'[r],r]+(((r^3)(Mv^2))/(((((lf1[r])^2)+(r^2)))^2))V[r]==0;

solv[virval_,mvval_]:=Module[{rule},
rule={vir->virval,Mv->mvval};NDSolve[{Vode/.rule,V[vir/.rule]==1,
V'[vir/.rule]==0},V,{r,vir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;

Vuv[vir_?NumericQ,Mv_?NumericQ]:=(V[rmax]/.solv[vir,Mv]  );
Vmass1 = i/.FindRoot[Vuv[mir1,i],{i,mir1,0,rmax}] ;

Print["(Subscript[m, V])=",Vmass1];

KVode:=D[(r^3)KV'[r],r]==0;
solkv[kvirval_]:=Module[{rule},
rule={kvir->kvirval};NDSolve[{KVode/.rule,KV[kvir/.rule]==1,
KV'[kvir/.rule]==0},KV,{r,kvir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;

NormV := Evaluate[(Sqrt[NIntegrate[ (r^3)(((V[r]/.solv[mir1,Vmass1]))^2)/((((r^2)+((lf1[r])^2)))^2)/g5sq, {r,mir1,rmax}]])]

vode[r_]:=Evaluate[(V[r]/.solv[mir1,Vmass1])/ NormV];
vodep[r_]:=vode'[r];
vodepp[r_]:=vode''[r];

kvode[r_]:=Evaluate[(KV[r]/.solkv[mir1])(Sqrt[((g5sq* Nfp *dr2)/(24 (Pi^2)))])/(KV[If[100*mir1>rmax,rmax,100*mir1]]/.solkv[mir1])];
fvsq1=(1/g5sq)NIntegrate[(-3 (r^2) vodep[r]-(r^3) vodepp[r])kvode[r],{r,mir1,If[100*mir1>rmax,rmax,100*mir1]}];
Print["(Subscript[f, V])=",(Sqrt[fvsq1])];

piode:=D[(r^3)lf1[r]^2Π'[r],r]+(((r^3)lf1[r]^2 (Mv^2))/(((((lf1[r])^2)+(r^2)))^2))Π[r]==0;

solpi[virval_,mvval_]:=Module[{rule},
rule={vir->virval,Mv->mvval};NDSolve[{piode/.rule,Π[vir/.rule]==1,
Π'[vir/.rule]==0},Π,{r,vir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;

piuv[vir_?NumericQ,Mv_?NumericQ]:=(Π[rmax]/.solpi[vir,Mv]  );
pimass1 = i/.FindRoot[piuv[mir1,i],{i,mir1,0,rmax}] ;
Print["(Subscript[m, π])=",pimass1];

pimass1=0.0

ode2:=D[(r^3)Π'[r],r]-g5sq ((((lf1[r])^2)(r^3)Π[r])/((((r^2)+((lf1[r])^2)))^2))==0;
solp[ival_]:=Module[{rule},rule={k->ival};NDSolve[{ode2/.rule, Π[k/.rule]==mir1,
Π'[k/.rule]==0},Π,{r,k/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First];
pode[r_]:=Evaluate[(Π[r]/.solp[mir1])(Sqrt[((g5sq* Nfp* dr2)/(24 (Pi^2)))])/(Π[If[100*mir1>rmax,rmax,100*mir1]]/.solp[mir1])];
podep[r_]:=pode'[r];
podepp[r_]:=pode''[r];
fpisq1=-(1/g5sq)NIntegrate[(-3 (r^2) podep[r]-(r^3) podepp[r])pode[r],{r,mir1,If[100*mir1>rmax,rmax,100*mir1]}];

Print["(Subscript[f, π])=",(Sqrt[fpisq1])];

Smax[Ms1_?NumericQ]:=Module[{sltnS},sltnS=NDSolve[{D[(r^3)S'[r],r]+(((r^3)(Ms1^2))/(((((lf1[r])^2)+(r^2)))^2))S[r]-del1[Log[(Sqrt[(r^2)+((lf1[r])^2)])]]S[r]r-((r ((lf1[r])^2))/((((lf1[r])^2)+(r^2))))del1'[Log[(Sqrt[(r^2)+((lf1[r])^2)])]]S[r]==0,S[mir1]==1,S'[mir1]==0},S,{r,mir1,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13];(S[r]/.sltnS/.r->rmax)];

Ms1Sol=FindRoot[Smax[Ms1]==0,{Ms1,0}]; (*adjust to avoid oscillatory behaviour!*)
Ms1=Ms1/.Ms1Sol;

numpts=10;
Sode:=D[(r^3)S'[r],r]+(((r^3)(Ms1^2))/(((((lf1[r])^2)+(r^2)))^2))S[r]-del1[Log[(Sqrt[(r^2)+((lf1[r])^2)])]]S[r]r-((r ((lf1[r])^2))/((((lf1[r])^2)+(r^2))))del1'[Log[(Sqrt[(r^2)+((lf1[r])^2)])]]S[r]==0;

sltnS=NDSolve[{Sode,S[mir1]==1,S'[mir1]==0},S,{r,mir1,rmax}];
Plot[S[r]/.sltnS,{r,mir1,0.00001rmax} ,BaseStyle->{18,FontFamily->"Times New Roman"},AxesLabel->{"ρ","S(ρ)"},PlotRange->{{mir1,0.00001rmax},{-0.3,1.1}}, PlotStyle->{Thick,Orange}]
S[r]/.sltnS/.r:>rmax
Smass1=Ms1;
Print["(Subscript[m, σ])=",Smass1];

KSode:=D[(r^3)KS'[r],r]-del1[Log[(Sqrt[(r^2)+(((lf1[r]))^2)])]]KS[r]r-((r ((lf1[r])^2))/((((lf1[r])^2)+(r^2))))del1'[Log[(Sqrt[(r^2)+(((lf1[r]))^2)])]]KS[r]==0;
sols[virval_]:=Module[{rule},
rule={vir->virval};NDSolve[{Sode/.rule,S[vir/.rule]==1,
S'[vir/.rule]==0},S,{r,vir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;
solkv[kvirval_]:=Module[{rule},
rule={kvir->kvirval};NDSolve[{KSode/.rule,KS[kvir/.rule]==1,
KS'[kvir/.rule]==0},KS,{r,kvir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;
NormS := Evaluate[(Sqrt[NIntegrate[ r^3 ((((S[r]/.sols[mir1]))^2)/(((r^2+((lf1[r])^2)))^2)), {r,mir1,rmax}]])]
Soode[r_]:=Evaluate[(S[r]/.sols[mir1])/ NormS];
(*Soode[r_]:=Evaluate[(S[r]/.sltnS)/ NormS];*)
Soodep[r_]:=Soode'[r];
Soodepp[r_]:=Soode''[r];
ksode[r_]:=Evaluate[(KS[r]/.solkv[mir1])(Sqrt[((Nfp*dr2)/(24 (Pi^2)))])/(KS[If[100*mir1>rmax,rmax,100*mir1]]/.solkv[mir1])];
fssq1=NIntegrate[(-3 (r^2) Soodep[r]-(r^3) Soodepp[r])ksode[r],{r,mir1,If[100*mir1>rmax,rmax,100*mir1]}];
Print["(Subscript[f, σ])=",(Sqrt[fssq1])];

Amax[Ma1_?NumericQ]:=Module[{sltnA},sltnA=NDSolve[{D[(r^3)A'[r],r]+(((r^3)((Ma1^2)-g5sq ((lf1[r])^2)))/((((((lf1[r])^2)+(r^2)))^2) ))A[r]==0,A[mir1]==1,A'[mir1]==0},A,{r,mir1,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13];
(A[r]/. sltnA)/. r->rmax];

Ma1Sol=FindRoot[Amax[Ma1]==0,{Ma1,3.6*^-7}]; (*adjust to avoid oscillatory behaviour!*)
Ma1=Ma1/. Ma1Sol

Aode:=D[(r^3)A'[r],r]+(((r^3)((Ma1^2)-g5sq ((lf1[r])^2)))/(((((lf1[r])^2)+(r^2)))^2))A[r]==0;
sltnA=NDSolve[{Aode,A[mir1]==1,A'[mir1]==0},A,{r,mir1,rmax}];
Plot[A[r]/.sltnA,{r,mir1,0.000001rmax} ,BaseStyle->{18,FontFamily->"Times New Roman"},AxesLabel->{"ρ","A(ρ)"},PlotRange->{{mir1,0.000001rmax},{-0.3,1.1}}, PlotStyle->{Thick,Orange}]
table1=Table[A[ρ]/.sltnA,{ρ,mir1,rmax,(10^(-1))}]//Flatten;
AllTrue[table1,#≥0&]
A[r]/.sltnA/.r:>rmax
Amass1=Ma1;
Print["(Subscript[m, A])=",Amass1];

Athing[r_] := First[Evaluate[A[r] /. sltnA]]

NormA := Evaluate[(Sqrt[NIntegrate[ (r^3)((Athing[r])^2)/((((r^2)+((lf1[r])^2)))^2)/g5sq, {r,mir1,rmax}]])]

apode[r_]:=Evaluate[(Athing[r])/NormA];
apodep[r_]:=apode'[r];
apodepp[r_]:=apode''[r];
fasq1=(1/g5sq)NIntegrate[(-3 (r^2) apodep[r]-(r^3) apodepp[r])pode[r],{r,mir1,If[100*mir1>rmax,rmax,100*mir1]}];

Print["(Subscript[f, A])=",(Sqrt[fasq1])];

Fundamental

Nfp=nfr1;(*number of flavours used to calculate fπ*)
g5sq:=((12 Pi^2)/(dr1 nfr1));
α2:=(π/(3.*c2r2)); (*Value of α[1] compare to output at end to tune for δ*)
B0:= (1/(6 π))(11 c2g -4tr1 nfr1 -4 tr2 nfr2) ;(*one-loop*) (*doubled factors to 4 because Dirac fermions*)
B1:=(1/(24 (π^2)))(34(c2g^2)- (20c2g + 12c2r1)tr1 nfr1- (20c2g + 12c2r2)tr2 nfr2 ) ;(*two-loop*) (*we have doubled the factors to 20 and 12*)
eqalp2:=alp2'[mu2] + (B0 alp2[mu2]^2  + B1 alp2[mu2]^3)==0
sol2 := NDSolve [{ eqalp2, alp2[0]==α2}, alp2, {mu2,-100, (10^24)},Method->{"StiffnessSwitching"}];
alph2[mu2_]:= First[Evaluate[alp2[mu2]/. sol2]]
αv2[mu_]:=alph2[mu]
γ2[mu2_]:=(3/(2π)) *c2r1*αv2[mu2] 
δ2[mu2_]:=-2γ2[mu2]
δ2[0]
del2[mu2_]=-2γ2[mu2];

plotdel2v2=Plot[del2[μ+Log[mρ]],{μ,-100,100},PlotStyle->{Thickness[0.009],Red},Frame->True,FrameLabel->{Style["ln (μ/(Subscript[m, ρ, 15]))",Black,17],Style["((Δm)^2)",Black,17]},FrameTicksStyle->Directive[Black,17],PlotTheme->"Scientific",FrameStyle->{{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black}},BaseStyle->{17,FontFamily->"Times New Roman"},AxesStyle->{Thick,Black},PerformanceGoal->"Quality",PlotRange->{{-60,60},{-4,0}}]

muv=0.0mρ;(*UV mass (IF using IR mass comment this and the line to calculate mir later)*)
guess2=1.17*(10^(-13));
ode:=D[(r^3)L'[r],r]-δ2[Log[Sqrt[(r^2)+((L[r])^2)] ]]r L[r]==0
solm[v_]:= NDSolve[{ode,L[v ]==v,
L'[v]==0},L,{r,v,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First
LinfBm[ v_?NumericQ ]:=(L[rmax]/.solm[v])-muv
(*Find mir from muv (comment out if using mir*)
mir2= m/.FindRoot[LinfBm[m],{m,guess2}]
lf2[ρ_]:=Evaluate[L[ρ]/.solm[ mir2]]
plot2=Show[Plot[((lf2[r*mρ])/(mρ)),{r,(mir2/(mρ)),rmax},PlotStyle->{Thickness[0.009],Red},Frame->True,FrameLabel->{Style["(ρ/(Subscript[m, ρ, 15]))",Black,Bold,17],Style["((L ((Subscript[ρm, ρ, 15])))/(Subscript[m, ρ, 15]))",Black,Bold,17]},FrameTicksStyle->Directive[Black,17],PlotTheme->"Scientific",FrameStyle->{{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black}},BaseStyle->{17,FontFamily->"Times New Roman"},AxesStyle->{Thick,Black},PerformanceGoal->"Quality", PlotRange->{{0,rmax},{-0.01,0.4}},PlotLegends->"Red: Vacuum for the fundamental"],Plot[x,{x,0,5}, PlotStyle->{Thickness[0.009],Black}]]
(*The above is the "embedding"*)
(*To get the point that we set the wall*)
(*It should be the IR-mass, but check for better or for worse*)
{r,lf2[r]}/.FindRoot[lf2[r]==r,{r,guess2}]
(*Now, we can write L(ρ=(Subscript[ρ, IR]))=(Subscript[ρ, IR]) and L'(ρ=(Subscript[ρ, IR]))=0*)
table1=Table[lf2[r],{r,mir2,rmax,(10^(-2))}];
AllTrue[table1,#>0&]
lf1[r]/.r:>10//SetPrecision[#,2]&
L2[ρ_]=lf2[ρ];

SU(2) with only the fundamental below log((Sqrt[2])(Subscript[m, IR,G]))

Log[Sqrt[2]*mir1];
αv2[Log[Sqrt[2]*mir1]]

Nfp=nfr1;(*number of flavours used to calculate fπ*)
g5sq:=((12 Pi^2)/(dr1 nfr1));
cut3:=Log[Sqrt[2]*mir1] (*cutoff value from log((Sqrt[2])(Subscript[m, IR]))*)
α3:=0.2187289174375807`; (*Value of α[1] at the log((Sqrt[2])(Subscript[m, IR])) cutoff compare to output at end to tune for δ*)
B0:=(1/(6π))(11 c2g -2 tr1 nfr1 );(*one-loop*)
B1:=(1/(24(π^2)))(34(c2g^2)- (20c2g + 12 c2r1)tr1 nfr1);(*two-loop*) (*we have doubled the factors to 20 and 12*)
eqalp3:=alp3'[mu3] + (B0 alp3[mu3]^2  + B1 alp3[mu3]^3)==0
sol3:= NDSolve [{ eqalp3, alp3[cut3]==α3}, alp3, {mu3,-100, cut3},Method->{"StiffnessSwitching"}];
alph3[mu3_]:= First[Evaluate[alp3[mu3]/. sol3]]
αv3[mu_]:=alph3[mu]
γ3[mu3_]:=(3/(2π)) *c2r1*αv3[mu3] 
δ3[mu3_]:=-2γ3[mu3]
del3[mu3_]=-2γ3[mu3];

plotdel3v2=Plot[del3[x+Log[mρ]],{x,-80,cut3-Log[mρ]},PlotStyle->{Thickness[0.009],Green},Frame->True,FrameLabel->{Style["ln (μ/(Subscript[m, ρ, 15]))",Black,17],Style["((Δm)^2)",Black,17]},FrameTicksStyle->Directive[Black,17],PlotTheme->"Scientific",FrameStyle->{{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black}},BaseStyle->{17,FontFamily->"Times New Roman"},AxesStyle->{Thick,Black},PerformanceGoal->"Quality",PlotRange->{{-60,20},{-10.5,1}}]

Show the above three runnings

Show[plotdel1v1,plotdel2v2,plotdel3v2]

F____full model and then integrating out

del[x_]=Piecewise[{{del3[x],-100≤x≤cut3},{del2[x],cut3<x≤(10^24)}}]

Plot[del[x],{x,-100,100}]

Show[Plot[del[x],{x,-1.9,1}],Plot[del3[x],{x,-1.9,cut3},PlotStyle->Red],Plot[del2[x],{x,cut3,1},PlotStyle->Purple],PlotRange->All,ImageSize->Large]

muv=0.0mρ;(*UV mass (IF using IR mass comment this and the line to calculate mir later)*)
guess3=2.4*(10^(-8)); (*initial guess for IR mass*)
ode:=D[(r^3)L'[r],r]-del[Log[Sqrt[(r^2)+((L[r])^2)] ]]r L[r]==0
solm[v_]:= NDSolve[{ode,L[v ]==v,
L'[v]==0},L,{r,v,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity]//First
LinfBm[ v_?NumericQ ]:=(L[rmax]/.solm[v])-muv
(*Find mir from muv (comment out if using mir*)
Quiet[mirrandom = m/.FindRoot[LinfBm[m],{m,guess3},WorkingPrecision->100]]
Quiet[lfnew[ρ_]:=Evaluate[L[ρ]/.solm[ mirrandom]]]
plot3=Show[Plot[((lfnew[r*mρ])/(mρ)),{r,(mirrandom/(mρ)),rmax},PlotStyle->{Thick,Orange},BaseStyle->{18,FontFamily->"Times New Roman"},AxesLabel->{"ρ","L(ρ)"}, PlotRange->{{0,rmax},{0,0.4}},PlotLegends->"Orange: Vacuum for the 2-index"],Plot[x,{x,0,0.4}, PlotStyle->{Thick,Black}]]

(*The above is the "embedding"*)
(*To get the point that we set the wall*)
(*It should be the IR-mass, but check for better or for worse*)
{r,lfnew[r]}/.FindRoot[lfnew[r]==r,{r,guess3}]
(*Now, we can write L(ρ=(Subscript[ρ, IR]))=(Subscript[ρ, IR]) and L'(ρ=(Subscript[ρ, IR]))=0*)

ode:=D[(r^3)L'[r],r]-del[Log[Sqrt[(r^2)+((L[r])^2)] ]]r L[r]==0
v := mirrandom;
(*take from the IR mass of the piecewise vacuum*)
Quiet[sltn=NDSolve[{ode,L[v ]==v,
L'[v]==0},L,{r,v,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->7,AccuracyGoal->7]//First
]
randomplot3=Show[Plot[((L[ρ*mρ] /. sltn)/(mρ)),{ρ,(v/(mρ)),rmax},PlotStyle->{Thickness[0.009],Green},Frame->True,FrameLabel->{Style["(ρ/(Subscript[m, ρ, 15]))",Black,Bold,17],Style["((L ((Subscript[ρm, ρ, 15])))/(Subscript[m, ρ, 15]))",Black,Bold,17]},FrameTicksStyle->Directive[Black,17],PlotTheme->"Scientific",FrameStyle->{{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black}},BaseStyle->{17,FontFamily->"Times New Roman"},AxesStyle->{Thick,Black},PerformanceGoal->"Quality", PlotRange->{{0,rmax},{-0.05,0.4}},PlotLegends->"Green: Piecewise vacuum"],Plot[x,{x,0,0.4}, PlotStyle->{Thickness[0.009],Black}]]
L[ρ]/.sltn/.ρ->10//SetPrecision[#,2]&
Linter[ρ_]=L[ρ]/.sltn

Show the above three embeddings

Show[plot1,plot2,randomplot3]

Define the vacuum for the equations of motion.

lf[ρ_]=Linter[ρ];
mir=v;
(*Plot[Linter[r], {r,v,10}, PlotRange -> {-1,1},PlotStyle->{Thick,Red}]
Plot[lf[r], {r,v,10}, PlotRange -> {-1,1},PlotStyle->{Thick,Red}]*)
mir

The equations of motion - Fundamentals

(* The automatic code works here *) 
Vode:=D[(r^3)V'[r],r]+(((r^3)(Mv^2))/(((((lf[r])^2)+(r^2)))^2))V[r]==0;

solv[virval_,mvval_]:=Module[{rule},
rule={vir->virval,Mv->mvval};NDSolve[{Vode/.rule,V[vir/.rule]==1,
V'[vir/.rule]==0},V,{r,vir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;

Vuv[vir_?NumericQ,Mv_?NumericQ]:=(V[rmax]/.solv[vir,Mv]  );
Vmass = i/.FindRoot[Vuv[mir,i],{i,mir,0,rmax}] ;

Print["(Subscript[m, V])=",Vmass];

(Vmass1/Vmass)

KVode:=D[(r^3)KV'[r],r]==0;
solkv[kvirval_]:=Module[{rule},
rule={kvir->kvirval};NDSolve[{KVode/.rule,KV[kvir/.rule]==1,
KV'[kvir/.rule]==0},KV,{r,kvir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;

NormV := Evaluate[(Sqrt[NIntegrate[ (r^3)(((V[r]/.solv[mir,Vmass]))^2)/((((r^2)+((lf[r])^2)))^2)/g5sq, {r,mir,rmax}]])]

vode[r_]:=Evaluate[(V[r]/.solv[mir,Vmass])/ NormV];
vodep[r_]:=vode'[r];
vodepp[r_]:=vode''[r];

kvode[r_]:=Evaluate[(KV[r]/.solkv[mir])(Sqrt[((g5sq* Nfp *dr1)/(24 (Pi^2)))])/(KV[If[100*mir>rmax,rmax,100*mir]]/.solkv[mir])];
fvsq=(1/g5sq)NIntegrate[(-3 (r^2) vodep[r]-(r^3) vodepp[r])kvode[r],{r,mir,If[100*mir>rmax,rmax,100*mir]}];
Print["(Subscript[f, V])=",(Sqrt[fvsq])];

piode:=D[(r^3)((lf[r])^2)Π'[r],r]+(((r^3)((lf[r])^2)(Mv^2))/(((((lf[r])^2)+(r^2)))^2))Π[r]==0;

solpi[virval_,mvval_]:=Module[{rule},
rule={vir->virval,Mv->mvval};NDSolve[{piode/.rule,Π[vir/.rule]==1,
Π'[vir/.rule]==0},Π,{r,vir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;

piuv[vir_?NumericQ,Mv_?NumericQ]:=(Π[rmax]/.solpi[vir,Mv]  );
pimass = i/.FindRoot[piuv[mir,i],{i,mir,0,rmax}] ;
Print["(Subscript[m, π])=",pimass];

((((pimass/(mρ))))^2)

ode2:=D[(r^3)Π'[r],r]-g5sq ((((lf[r])^2)(r^3)Π[r])/((((r^2)+((lf[r])^2)))^2))==0;
solp[ival_]:=Module[{rule},rule={k->ival};NDSolve[{ode2/.rule, Π[k/.rule]==mir,
Π'[k/.rule]==0},Π,{r,k/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First];
pode[r_]:=Evaluate[(Π[r]/.solp[mir])(Sqrt[((g5sq* Nfp* dr1)/(24 (Pi^2)))])/(Π[If[100*mir>rmax,rmax,100*mir]]/.solp[mir])];
podep[r_]:=pode'[r];
podepp[r_]:=pode''[r];
fpisq=-(1/g5sq)NIntegrate[(-3 (r^2) podep[r]-(r^3) podepp[r])pode[r],{r,mir,If[100*mir>rmax,rmax,100*mir]}];

Print["(Subscript[f, π])=",(Sqrt[fpisq])];

mirf=mir+0.*(10^(-8)); (*adjust to avoid singularity when root-finding*)

Smax[Ms_?NumericQ]:=Module[{sltnS},sltnS=NDSolve[{D[(r^3)S'[r],r]+(((r^3)(Ms^2))/(((((lf[r])^2)+(r^2)))^2))S[r]-del[Log[(Sqrt[(r^2)+((lf[r])^2)])]]S[r]r-((r ((lf[r])^2))/((((lf[r])^2)+(r^2))))del'[Log[(Sqrt[(r^2)+((lf[r])^2)])]]S[r]==0,S[mirf]==1,S'[mirf]==0},S,{r,mirf,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13];(S[r]/.sltnS/.r->rmax)];

MsSol=FindRoot[Smax[Ms]==0,{Ms,8.6*^-9}]; (*adjust to avoid oscillatory behaviour!*)
Ms=Ms/.MsSol;

numpts=10;

Sode:=D[(r^3)S'[r],r]+(((r^3)(Ms^2))/(((((lf[r])^2)+(r^2)))^2))S[r]-del[Log[(Sqrt[(r^2)+(((lf[r]))^2)])]]S[r]r-((r ((lf[r])^2))/((((lf[r])^2)+(r^2))))del'[Log[Sqrt[(r^2)+(((lf[r]))^2)]]]S[r]==0;

Quiet[sltnS=NDSolve[{Sode,S[mirf]==1,S'[mirf]==0},S,{r,mirf,rmax}]];
Plot[S[r]/.sltnS,{r,mirf,0.000001rmax} ,BaseStyle->{18,FontFamily->"Times New Roman"},AxesLabel->{"ρ","S(ρ)"},PlotRange->{{mirf,0.000001rmax},{-0.3,1.1}}, PlotStyle->{Thick,Orange}]
S[r]/.sltnS/.r:>rmax
Smass=Ms;
Print["(Subscript[m, σ])=",Smass];

mirf=mir+0.*(10^(-8));
KSode:=D[(r^3)KS'[r],r]-del[Log[Sqrt[(r^2)+(((lf[r]))^2)]]]KS[r]r-((r ((lf[r])^2))/((((lf[r])^2)+(r^2))))del'[Log[Sqrt[(r^2)+(((lf[r]))^2)]]]KS[r]==0;
sols[virval_]:=Module[{rule},
rule={vir->virval};NDSolve[{Sode/.rule,S[vir/.rule]==1,
S'[vir/.rule]==0},S,{r,vir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;
solkv[kvirval_]:=Module[{rule},
rule={kvir->kvirval};NDSolve[{KSode/.rule,KS[kvir/.rule]==1,
KS'[kvir/.rule]==0},KS,{r,kvir/.rule,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13]//First] ;
Quiet[NormS := Evaluate[Sqrt[NIntegrate[ r^3  (S[r]/.sols[mirf])^2/(r^2+lf[r]^2)^2, {r,mirf,rmax}]]]]
Quiet[Soode[r_]:=Evaluate[(S[r]/.sols[mirf])/ NormS]];
Soodep[r_]:=Soode'[r];
Soodepp[r_]:=Soode''[r];
Quiet[ksode[r_]:=Evaluate[(KS[r]/.solkv[mirf])*(Sqrt[((Nfp*dr1)/(24 (Pi^2)))])/(KS[If[100*mirf>rmax,rmax,100*mirf]]/.solkv[mirf])]];
fssq=NIntegrate[(-3 (r^2) Soodep[r]-(r^3) Soodepp[r])ksode[r],{r,mirf,If[100*mirf>rmax,rmax,100*mirf]}];
Print["(Subscript[f, σ])=",(Sqrt[fssq])];

Amax[Ma_?NumericQ]:=Module[{sltnA},sltnA=NDSolve[{D[(r^3)A'[r],r]+(((r^3)((Ma^2)-g5sq ((lf[r])^2)))/((((((lf[r])^2)+(r^2)))^2) ))A[r]==0,A[mir]==1,A'[mir]==0},A,{r,mir,rmax},Method->{"StiffnessSwitching"},MaxSteps->Infinity,PrecisionGoal->13,AccuracyGoal->13];
(A[r]/. sltnA)/. r->rmax];

MaSolution=FindRoot[Amax[Ma]==0,{Ma,2.3*^-8}]; (*adjust to avoid oscillatory behaviour!*)
Ma=Ma/. MaSolution;

Aode:=D[(r^3)A'[r],r]+(((r^3)((Ma^2)-g5sq ((lf[r])^2)))/(((((lf[r])^2)+(r^2)))^2))A[r]==0;

sltnA=NDSolve[{Aode,A[mir]==1,A'[mir]==0},A,{r,mir,rmax}];
Plot[A[r]/.sltnA,{r,mir,0.000001rmax} ,BaseStyle->{18,FontFamily->"Times New Roman"},AxesLabel->{"ρ","A(ρ)"},PlotRange->{{mir,0.000001rmax},{-0.1,1.1}}, PlotStyle->{Thick,Orange}]
table1=Table[A[ρ]/.sltnA,{ρ,mir,rmax,(10^(-1))}]//Flatten;
AllTrue[table1,#≥0&]
A[r]/.sltnA/.r:>rmax
Amass=Ma;
Print["(Subscript[m, A])=",Amass];

Athing[r_] := First[Evaluate[A[r] /. sltnA]]

NormA := Evaluate[Sqrt[NIntegrate[ r^3  Athing[r]^2/(r^2+lf[r]^2)^2/g5sq, {r,mir,rmax}]]]

apode[r_]:=Evaluate[(Athing[r])/NormA];
apodep[r_]:=apode'[r];
apodepp[r_]:=apode''[r];
fasq=(1/g5sq)NIntegrate[(-3 (r^2) apodep[r]-(r^3) apodepp[r])pode[r],{r,mir,If[100*mir>rmax,rmax,100*mir]}];

Print["(Subscript[f, A])=",(Sqrt[fasq])];

Results as a table.

Print["Rescaled"];
Print["(Subscript[m, ρ, S])=",(Vmass1/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[m, ρ, F])=",(Vmass/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[m, σ, S])=",(Smass1/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[m, σ, F])=",(Smass/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[m, A, S])=",(Amass1/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[m, A, F])=",(Amass/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, ρ, S])=",((Sqrt[fvsq1])/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, ρ, F])=",((Sqrt[fvsq])/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, π, S])=",((Sqrt[fpisq1])/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, π, F])=",((Sqrt[fpisq])/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, σ, S])=",((Sqrt[fssq1])/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, σ, F])=",((Sqrt[fssq])/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, A, S])=",((Sqrt[fasq1])/Vmass1)//SetPrecision[#,3]&];
Print["(Subscript[f, A, F])=",((Sqrt[fasq])/Vmass1)//SetPrecision[#,3]&];

Print["Unscaled"];
Print["(Subscript[m, ρ, S])=",Vmass1//SetPrecision[#,3]&];
Print["(Subscript[m, ρ, F])=",Vmass//SetPrecision[#,3]&];
Print["(Subscript[m, σ, S])=",Smass1//SetPrecision[#,3]&];
Print["(Subscript[m, σ, F])=",Smass//SetPrecision[#,3]&];
Print["(Subscript[m, A, S])=",Amass1//SetPrecision[#,3]&];
Print["(Subscript[m, A, F])=",Amass//SetPrecision[#,3]&];
Print["(Subscript[f, ρ, S])=",(Sqrt[fvsq1])//SetPrecision[#,3]&];
Print["(Subscript[f, ρ, F])=",(Sqrt[fvsq])//SetPrecision[#,3]&];
Print["(Subscript[f, π, S])=",(Sqrt[fpisq1])//SetPrecision[#,3]&];
Print["(Subscript[f, π, F])=",(Sqrt[fpisq])//SetPrecision[#,3]&];
Print["(Subscript[f, σ, S])=",(Sqrt[fssq1])//SetPrecision[#,3]&];
Print["(Subscript[f, σ, F])=",(Sqrt[fssq])//SetPrecision[#,3]&];
Print["(Subscript[f, A, S])=",(Sqrt[fasq1])//SetPrecision[#,3]&];
Print["(Subscript[f, A, F])=",(Sqrt[fasq])//SetPrecision[#,3]&];