SetOptions[LogLinearPlot,AxesOrigin->Automatic,BaseStyle->{FontSize->18},FrameTicksStyle->Directive[16],Frame->True,FrameStyle->Darker[Gray,1/2],LabelStyle->{Italic,FontSize->20}];
DarkGreen=Darker[Green,1/3];

Input parameters

(* beta function coefficients *)
betas={CF->4/3,CA->3,nf->5,b0->11-2/3*nf,b1->102-38/3nf,b2->2857/54*CA^3+(CF^2-205/18 CF*CA-1415/54*CA^2) nf+(11/9*CF+79/54*CA)*nf^2,b3->149753/6+3564 Zeta[3]+(-1078361/162-6508/27 Zeta[3]) nf+(50065/162+6472/81 Zeta[3]) nf^2+1093/729 nf^3};

(* Initial condition for alphas *)
lowmu0={mu0->91187/1000,asmu0->118/1000};

Strong Coupling Evolution (numerical)

(* Routine which solves the differential equation for alpha_s with the beta-function known at n loops *)

numsol[n_]:=NDSolve[{(mu asnum'[mu]+asnum[mu]^2/(2Pi)  (Sum[b[i] (asnum[mu]/4/Pi)^i,{i,0,n}])==0)/.b[0]->b0/.b[1]->b1/.b[2]->b2/.b[3]->b3//.betas,(asnum[mu0]==asmu0)/.lowmu0},asnum,{mu,1/2,100}]

asnum0=asnum[mu]/.numsol[0];
asnum1=asnum[mu]/.numsol[1];
asnum2=asnum[mu]/.numsol[2];
asnum3=asnum[mu]/.numsol[3];

asnumeric=LogLinearPlot[{asnum0,asnum1,asnum2,asnum3},{mu,1/2,100},PlotStyle->{{Dashed,Black},{Dashed,Blue},{Dashed,Red},{Dashed,DarkGreen}},Frame->True,FrameLabel->{"μ (GeV)","(Subscript[α, s])"},ImageSize->Large]

Strong Coupling Evolution (analytical)

(* the analytic solution of the equation for alpha_s can be obtained perturbatively. We first expand the integrand in alpha_s, by means of the variable t=-2Pi/beta0/alpha_s *)
 
den[t_]=1-b1/2/b0^2 1/t+b2/4/b0^3/t^2-b3/8/b0^4/t^3;

(* Compute the bihat coefficients of function G(t) by expanding the integrand denominator in 1/t (see the definition of the bihat coefficients in the exercise notes) *)
b1hat=SeriesCoefficient[Series[1/den[t],{t,Infinity,3}],1]
b2hat=SeriesCoefficient[Series[1/den[t],{t,Infinity,3}],2]
b3hat=SeriesCoefficient[Series[1/den[t],{t,Infinity,3}],3]

(* Equation for x=t0/t to be solved iteratively *)
xeq[x_]:= Xf + eps b1hat/t0 (-Log[x])-eps^2 b2hat/t0^2(1-1/x)-eps^3 b3hat/2 /t0^3(1-1/x^2)
(* The term Xf is O(1), and terms proportional to 1/t0 are O(alpha_s). To keep track of the perturbative order I introduce an auxiliary "small" parameter eps. Make eps=1 after the expansion *)

LL solution

(* The LL solution, x[0], is obtaineb by keeping the leading term in the equation xeq[x] *)

x[0]=Series[xeq[x],{eps,0,0}]//Normal

NLL solution

(* The NLL solution, x[1], is obtained by keeping up to NLL terms in the equation xeq evaluated at x[0]  *)

x[1]=Series[xeq[x[0]],{eps,0,1}]//Normal
delta[1]=Simplify[x[1]-x[0]]/.eps->1

NNLL solution

(* The NNLL solution, x[2], is obtained by keeping up to NNLL terms in the equation xeq evaluated at x[1]  *)

x[2]=Series[xeq[x[1]],{eps,0,2}]//Normal
delta[2]=Simplify[x[2]-x[1]]/.eps->1

N3LL solution

x[3]=Series[xeq[x[2]],{eps,0,3}]//Normal
delta[3]=Simplify[x[3]-x[2]]/.eps->1

(* Analytic solution for alpha_s(mu)=alpha_s(mu0)/x[i] at order i *)

as[i_,mu_]:=(asmu0/x[i])/.t0->-2Pi/b0/asmu0/.Xf->1+asmu0 b0/(2Pi)Log[mu/mu0]/.eps->1//.betas

asanalytic=LogLinearPlot[{as[0,mu]/.lowmu0,as[1,mu]/.lowmu0,as[2,mu]/.lowmu0,as[3,mu]/.lowmu0},{mu,1/2,100},PlotStyle->{Black,Blue,Red,DarkGreen},Frame->True,FrameLabel->{"μ (GeV)","(Subscript[α, s])"},ImageSize->Large]

Analytic (solid lines) vs numeric (dashed lines) results

Show[asnumeric,asanalytic,Frame->True,FrameLabel->{"μ (GeV)","(Subscript[α, s])"},ImageSize->Large,PlotRange->{{Log[1/2],Log[10]},{0.15,0.8}}]

diff0[mu_]:=(asnum0-(as[0,mu]/.lowmu0))/(as[0,mu]/.lowmu0)100
diff1[mu_]:=(asnum1-(as[1,mu]/.lowmu0))/(as[1,mu]/.lowmu0)100
diff2[mu_]:=(asnum2-(as[2,mu]/.lowmu0))/(as[2,mu]/.lowmu0)100
diff3[mu_]:=(asnum3-(as[3,mu]/.lowmu0))/(as[3,mu]/.lowmu0)100

LogLinearPlot[{diff0[mu],diff1[mu],diff2[mu],diff3[mu]},{mu,1/2,10},PlotStyle->{Black,Blue,Red,DarkGreen},Frame->True,FrameLabel->{"μ (GeV)","[((Subscript[α, s]))^num-((Subscript[α, s]))^ana]/((Subscript[α, s]))^ana (%)"},ImageSize->Large,PlotRange->{{Log[1/2],Log[10]},{0.0001,4}}]

Lambda QCD

G[t_]:=t + eps b1hat  Log[-t]-eps^2 b2hat/t-eps^3 b3hat/2 /t^2

(* Routine which computes LambdaQCD at i-th order, with the analytic result for alpha_s(mu) obtained at order asorder *) 

LambdaQCD[i_,mu_,asorder_]:=mu (Series[Exp[G[t]],{eps,0,i}]//Normal)/.eps->1/.t->-2Pi/b0/as[asorder,mu]//.betas

(* mu-dependence of LambdaQCD at at i-th order, with alpha_s evaluated at the same order *)

LogLinearPlot[{LambdaQCD[0,mu,0]/.lowmu0,LambdaQCD[1,mu,1]/.lowmu0,LambdaQCD[2,mu,2]/.lowmu0,LambdaQCD[3,mu,3]/.lowmu0},{mu,1/2,100},PlotStyle->{Black,Blue,Red,DarkGreen},FrameLabel->{"μ (GeV)","(Subscript[Λ, QCD]) (GeV)"},ImageSize->Large]

(* mu-dependence of LambdaQCD at at given order, with alpha_s evaluated at the 4th order (4-loop) *)

LogLinearPlot[{LambdaQCD[0,mu,3]/.lowmu0,LambdaQCD[1,mu,3]/.lowmu0,LambdaQCD[2,mu,3]/.lowmu0,LambdaQCD[3,mu,3]/.lowmu0},{mu,1/2,100},PlotStyle->{Black,Blue,Red,DarkGreen},FrameLabel->{"μ (GeV)","(Subscript[Λ, QCD]) (GeV)"},ImageSize->Large]