ρ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*^-15},
a[t],{t,0,500 Gyr}, Method->{"ImplicitRungeKutta","DifferenceOrder"->20}, WorkingPrecision->32]

α[t_]:=Evaluate[a[t]/.sol]

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[α[tt]-1,{tt,10*^9 yr}]; {"t0"->t0/Gyr "Gyr"}
N[{"ΩR"->ΩR,"ΩM"->ΩM,"ΩK"->ΩK,"ΩΛ"->ΩΛ}]
N[{"ρc"->ρc"kg/m³"}]

plot1=Plot[{α[t yr],(Sqrt[2]) ((ΩR)^(1/4))(Sqrt[H0 t yr]),1.01 ((((3/2)))^(2/3)) ((ΩM)^(1/3))(((H0 t yr))^(2/3)),(E^(H0 Sqrt[ΩΛ] (t-12.7*^9) yr))}, 
{t,0,40*^9}, GridLines->{{t0/yr}, {1}}, Frame->True, PlotRange->{All,All}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta,th}, (* Ω *)
{Red,th, Dashed}, (* Ωr *) 
{Blue,th,Dashed}, (* Ωm *)
{Darker@Green,th,Dashed}(* ΩΛ *)
},FrameTicksStyle->st]

plot2=LogLogPlot[{α[t yr],(Sqrt[2]) ((ΩR)^(1/4))(Sqrt[H0 t yr]),1.01 ((((3/2)))^(2/3)) ((ΩM)^(1/3))(((H0 t yr))^(2/3)),(E^(H0 Sqrt[ΩΛ] (t-12.7*^9) yr))}, 
{t, 1*^-10,500*^9}, GridLines->{{t0/yr}, {1}}, Frame->True, PlotRange->{All,All}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta,th}, (* Ω *)
{Red,th, Dashed}, (* Ωr *) 
{Blue,th,Dashed}, (* Ωm *)
{Darker@Green,th,Dashed}(* ΩΛ *)
},FrameTicksStyle->st]

plot3=LogPlot[{α[t yr],(Sqrt[2]) ((ΩR)^(1/4))(Sqrt[H0 t yr]),1.01 ((((3/2)))^(2/3)) ((ΩM)^(1/3))(((H0 t yr))^(2/3)),(E^(H0 Sqrt[ΩΛ] (t-12.7*^9) yr))},
 {t, 100000,40*^9}, GridLines->{{t0/yr}, {1}}, Frame->True, PlotRange->{All,All}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta,th}, (* Ω *)
{Red,th, Dashed}, (* Ωr *) 
{Blue,th,Dashed}, (* Ωm *)
{Darker@Green,th,Dashed}(* ΩΛ *)
},FrameTicksStyle->st]

plot4=LogLinearPlot[{α[t yr],(Sqrt[2]) ((ΩR)^(1/4))(Sqrt[H0 t yr]),1.01 ((((3/2)))^(2/3)) ((ΩM)^(1/3))(((H0 t yr))^(2/3)),(E^(H0 Sqrt[ΩΛ] (t-12.7*^9) yr))}, 
{t, 0,40*^9}, GridLines->{{t0/yr}, {1}}, Frame->True, PlotRange->{All,All}, 
ImagePadding->{{200, 2}, {60, 2}}, PlotPoints->im, ImageSize->im,
 AspectRatio->93/160, MaxRecursion->15, PlotStyle->
{{Darker@Magenta,th}, (* Ω *)
{Red,th, Dashed}, (* Ωr *) 
{Blue,th,Dashed}, (* Ωm *)
{Darker@Green,th,Dashed}(* ΩΛ *)
},FrameTicksStyle->st]

(* yukterez.net *)