ρr=((168132 )/100000)((8 (π^5) (kB^4) (T^4))/(15(c^5)(h^3))); (* Current Radiation+Neutrino Density *)
ρc=((3(H0^2))/(8π G)); (* Critical Density *)
H0=67150 ((m/sek)/Mpc); (* Current Hubble Parameter *)
T=(2725/1000) K; (* CMB Temperature *)
ΩR=((ρr)/(ρc)) ; (* Current Radiation+Neutrino Proportion *)

ΩM=315/1000;  (* Current Matter Proportion *)
ΩΛ=1-ΩR-ΩM;  (* Current Dark Energy Proportion *)
ΩK=1-ΩR-ΩM-ΩΛ;  (* Current Curvature Proportion *)

c=299792458m/sek;  (* Lightspeed *)
G=667384*^-16(m^3)(kg^(-1))(sek^(-2));  (* Gravitational Constant *)
kB=13806488*^-30 kg (m^2)/(sek^2)/K; (* Boltzmann Constant *)
h=662606957*^-42 kg (m^2)/sek; (* Planck Constant *)
Mpc=30856775777948584200000m;  (* Megaparsec *)
yr=36525*24*36*sek; (*Year *)
Gyr=10^9*yr; (* Gigayear *)
Glyr=Gyr*c; (* Gigalightyear *)
kg=m=sek=1; (* SI Units *)

sol=NDSolve[{((a'[t])/(a[t]))==H0 (Sqrt[((ΩR)/((a[t])^4))+((ΩM)/((a[t])^3))+((ΩK)/((a[t])^2))+ΩΛ]), a[0]==1*^-20},
a[t],{t,0,100 Gyr}, Method->{"ImplicitRungeKutta","DifferenceOrder"->20}, WorkingPrecision->32, MaxSteps->Infinity]
α[t_]:=Evaluate[(a[t]/.sol)[[1]]]

f=3 H0 t ((ΩM)^2);
a0[t_]:=If[t<τ1 yr, a1[t], If[t<τ2 yr, a2[t], a3[t]]];
a1[t_]:=(1/(4 ΩM))(2 (1-I (Sqrt[3])) ΩR+2 (1+I (Sqrt[3])) ΩR+(4 I (I+(Sqrt[3])) ((ΩR)^2))/((((f^2)-8f ((ΩR)^(3/2))+8 ((ΩR)^3)+(Sqrt[3]) ΩM (Sqrt[H0 t (f-8 ((ΩR)^(3/2)))]) Abs[f-4 ((ΩR)^(3/2))]))^(1/3))-I (-I+(Sqrt[3])) ((((f^2)-8f ((ΩR)^(3/2))+8 ((ΩR)^3)+(Sqrt[3]) ΩM (Sqrt[H0 t (f-8 ((ΩR)^(3/2)))]) Abs[f-4 ((ΩR)^(3/2))]))^(1/3)));
a2[t_]:=(1/(2 ΩM))(2 ΩR+(4 ((ΩR)^2))/((((f^2)-8f ((ΩR)^(3/2))+8 ((ΩR)^3)+(Sqrt[3]) ΩM (Sqrt[H0 t (f-8 ((ΩR)^(3/2)))]) Abs[f-4 ((ΩR)^(3/2))]))^(1/3))+((((f^2)-8f ((ΩR)^(3/2))+8 ((ΩR)^3)+(Sqrt[3]) ΩM (Sqrt[H0 t (f-8 ((ΩR)^(3/2)))]) Abs[f-4 ((ΩR)^(3/2))]))^(1/3)));
a3[t_]:=(((-1+(1/(ΩΛ))))^(1/3)) ((Sinh[(3/2) H0 t (Sqrt[ΩΛ])])^(2/3));

τ1=Chop[t/.FindRoot[a1[t]-a2[t], {t,1*^5 yr}, WorkingPrecision->32]]/yr
τ2=Chop[t/.FindRoot[a2[t]-a3[t], {t,1*^9 yr}, WorkingPrecision->32]]/yr

th= Thickness[0.007]; (* Plot Style *)
im=960; (* Image Width *)
st=Directive[Thick,FontFamily->"Lucida Console", FontSize->Round[48*(((im/1920))^(3/2))]];

t0=tt/.FindRoot[a0[tt]-1,{tt,10*^9 yr}]; {"t0"->t0/Gyr "Gyr"}
N[{"ΩR"->ΩR,"ΩM"->ΩM,"ΩK"->ΩK,"ΩΛ"->ΩΛ}]
N[{"ρc"->ρc"kg/m³"}]

t=τ yr;
plot[pf_,F_,τmin_,τmax_]:=pf[F, 
{τ,τmin,τmax}, GridLines->{{t0/yr, τ1, τ2}, {1}}, Frame->True, PlotRange->{All,All}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, WorkingPrecision->32, PlotStyle->
{{Darker@Magenta,th}, (* Ω *)
{Orange,th, Dashed}, (* Ωr *) 
{Red,th,Dashed}, (* Ωm *)
{Darker@Cyan,th,Dashed}(* ΩΛ *)
},FrameTicksStyle->st]

plot[Plot,{α[t],a1[t],a2[t],a3[t]},0,300000]

plot[Plot,{α[t],a1[t],a2[t],a3[t]},0,20*^9]

plot[LogLogPlot,{α[t],a1[t],a2[t],a3[t]},0.7576380963017885, 9.370690649675879*^10]

plot[LogLogPlot,{I,If[τ<τ1,a1[t],I],If[τ>τ1,If[τ<τ2,a2[t], I]],If[τ<τ2,I,a3[t]]},0.7576380963017885, 9.370690649675879*^10]

plot[LogLogPlot,{-1+α[t]/a0[t]},0.7576380963017885, 9.370690649675879*^10]