Ωr=((ρr)/(ρt)); (* Radiation+Neutrino Proportion *)

Ωm=((ρm)/(ρt)); (* Matter Proportion *)

Ωk=((ρk)/(ρt)); (* Curvature Proportion *)

ΩΛ=((ρΛ)/(ρt)); (* DE Proportion *)

ρr=((ρr0)/(a^4)); (* Radiation+Neutrino Density *)

ρm=((ρm0)/(a^3)); (* Matter Density *)

ρk=((ρk0)/(a^2)); (* Curvature Density *)

ρΛ=((ρΛ0)/(a^0)); (* DE Density *)

ρt=((3(H^2))/(8π G)); (* Critical Density *)

ρr0=((168132 )/100000)((8 (π^5) (kB^4) (T^4))/(15(c^5)(h^3))); (* Current Radiation+Neutrino Density *)

ρm0=Ωm0 ((3(H0^2))/(8π G)); (* Current Matter Density *)

ρk0=Ωk0 ((3(H0^2))/(8π G)); (* Current Curvature Density *)

ρΛ0=ΩΛ0 ((3(H0^2))/(8π G)); (* Current DE Density *)

H=H0(Sqrt[((Ωr0)/(a^4))+((Ωm0)/(a^3))+((Ωk0)/(a^2))+ΩΛ0]); (* Hubble Parameter *)

H0=67150 m/Mpc/sek; (* Current Hubble Parameter *)
T=2725/1000 K; (* CMB Temperature *)
Ωr0=ρr0 8π G/3/(H0^2); (* Current Radiation+Neutrino Proportion *)
Ωm0=315/1000;  (* Current Matter Proportion *)
ΩΛ0=1-Ωr0-Ωm0;  (* Current Dark Energy Proportion *)
Ωk0=1-Ωr0-Ωm0-ΩΛ0;  (* 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 *)
eV=1602176634*^-28kg (m^2)/(sek^2); (* Electronvolt *)
geV=1*^9 eV; (* Gigaelectronvolt *)
kg=m=sek=1; (* SI Units *)

sol=NDSolve[{((α'[t])/(α[t]))==H0 (Sqrt[((Ωr0)/((α[t])^4))+((Ωm0)/((α[t])^3))+((Ωk0)/((α[t])^2))+ΩΛ0]), α[0]==1*^-15},
α,{t,0,500 Gyr}, Method->{"ImplicitRungeKutta","DifferenceOrder"->20}, 
WorkingPrecision->32,MaxSteps->Infinity]

a:=Evaluate[(α[t]/.sol)[[1]]]
ȧ:=Evaluate[(α'[t]/.sol)[[1]]]
ä:=Evaluate[(α''[t]/.sol)[[1]]]
t0=t/.FindRoot[a-1,{t,10 Gyr}];
"t0"->t0/Gyr "Gyr"

th= Thickness[0.007]; (* Plot Style *)
im=960; (* Image Width *)
f=(a^3)/ρm0; (* Density to Energy Normalization *)
F=1; (* Plot Range Height Factor *)
t=time yr; (* Time Scale *)
st=Directive[Thick,FontFamily->"Lucida Console", FontSize->Round[48*(((im/1920))^(3/2))]];

o1=Plot[{Ωr,Ωm,ΩΛ}, {time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,All}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Ωr *) 
{Blue,th}, (* Ωm *)
{Darker@Green,th} (* ΩΛ *)
},FrameTicksStyle->st]
Export["o1.png", o1]

r1=Plot[{ρr (m^3)(c^2)/geV,ρm(m^3)(c^2)/geV,ρΛ(m^3)(c^2)/geV}, {time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True,PlotRange->{All,Automatic}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* ρr *) 
{Blue,th}, (* ρm *)
{Darker@Green,th} (* ρΛ *)
},FrameTicksStyle->st]
Export["r1.png", r1]

e1=Plot[{ρr f,ρm f,ρΛ f}, {time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True,PlotRange->{All,Automatic}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Er *) 
{Blue,th}, (* Em *)
{Darker@Green,th} (* EΛ *)
},FrameTicksStyle->st]
Export["e1.png", e1]

t1=Plot[{a}, {time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,Automatic} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta, th} (* a *)
},FrameTicksStyle->st]
Export["t1.png", t1]

o2=LogLogPlot[{Ωr,Ωm,ΩΛ}, {time, 1*^-6, 200*^9}, 
GridLines->{{t0/yr}, {0}}, Frame->True, PlotRange->{All,All} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Ωr *) 
{Blue,th}, (* Ωm *)
{Darker@Green,th} (* ΩΛ *)
},FrameTicksStyle->st]
Export["o2.png", o2]

r2=LogLogPlot[{ρr(m^3)(c^2)/geV,ρm(m^3)(c^2)/geV,ρΛ(m^3)(c^2)/geV}, {time, 1*^-6, 200*^9},
GridLines->{{t0/yr}, {0}}, Frame->True, PlotRange->{All,All} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* ρr *) 
{Blue,th}, (* ρm *)
{Darker@Green,th} (* ρΛ *)
},FrameTicksStyle->st]
Export["r2.png", r2]

e2=LogLogPlot[{ρr f,ρm f,ρΛ f}, {time, 1*^-6, 200*^9},
GridLines->{{t0/yr}, {0}}, Frame->True, PlotRange->{All,All} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Er *) 
{Blue,th}, (* Em *)
{Darker@Green,th} (* EΛ *)
},FrameTicksStyle->st]
Export["e2.png", e2]

t2=LogLogPlot[{a}, {time, 1*^-6, 200*^9},
GridLines->{{t0/yr}, {0}}, Frame->True, PlotRange->{All,All} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta, th} (* a *)
},FrameTicksStyle->st]
Export["t2.png", t2]

o3=LogPlot[{Ωr,Ωm,ΩΛ}, {time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,{1*^-5,2}} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Ωr *) 
{Blue,th}, (* Ωm *)
{Darker@Green,th} (* ΩΛ *)
},FrameTicksStyle->st]
Export["o3.png", o3]

r3=LogPlot[{ρr(m^3)(c^2)/geV,ρm(m^3)(c^2)/geV,ρΛ(m^3)(c^2)/geV},{time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,Automatic} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* ρr *) 
{Blue,th}, (* ρm *)
{Darker@Green,th} (* ρΛ *)
},FrameTicksStyle->st]
Export["r3.png", r3]

e3=LogPlot[{ρr f,ρm f,ρΛ f},{time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,Automatic} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Er *) 
{Blue,th}, (* Em *)
{Darker@Green,th} (* EΛ *)
},FrameTicksStyle->st]
Export["e3.png", e3]

t3=LogPlot[{a}, {time, 0, 25*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,Automatic} , MaxRecursion->15,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta, th} (* a *)
},FrameTicksStyle->st]
Export["t3.png", t3]

o4=LogLinearPlot[{Ωr,Ωm,ΩΛ},{time, 1, 1000*^9}, 
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,All} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Ωr *) 
{Blue,th}, (* Ωm *)
{Darker@Green,th} (* ΩΛ *)
},FrameTicksStyle->st]
Export["o4.png", o4]

r4=LogLinearPlot[{ρr(m^3)(c^2)/geV,ρm(m^3)(c^2)/geV,ρΛ(m^3)(c^2)/geV}, {time, 1*^8, 25*^9},
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,{-10/20,10}} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* ρr *) 
{Blue,th}, (* ρm *)
{Darker@Green,th} (* ρΛ *)
},FrameTicksStyle->st]
Export["r4.png", r4]

e4=LogLinearPlot[{ρr f,ρm f,ρΛ f}, {time, 1000, 25*^9},
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,{-10/20,10}} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Red,th}, (* Er *) 
{Blue,th}, (* Em *)
{Darker@Green,th} (* EΛ *)
},FrameTicksStyle->st]
Export["e4.png", e4]

t4=LogLinearPlot[{a},  {time, 1000, 25*^9},
GridLines->{{t0/yr}, {}}, Frame->True, PlotRange->{All,All} ,
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta, th} (* a *)
},FrameTicksStyle->st]
Export["t4.png", t4]