(* flat DE+matter universe*)

ΩM=315/1000;
ΩΛ=1-ΩM;
H0=217618*^-23;
Gyr=(10^7)*36525*24*3600;
Pc=30856775777948584;
a0=1;

(* scale factor *)

a[t_]:=(((-1+(1/(ΩΛ))))^(1/3)) ((Sinh[(3/2) H0 t (Sqrt[ΩΛ])])^(2/3));

(* t0 at a=1 *)

t0=NSolve[a[t ]==a0&&t>0,t][[1,1,2]]/Gyr; "t0"->t0 "Gyr"

(* scale factor plot *)

Plot[a[t Gyr],{t,0,2t0},Frame->True,GridLines->{{t0},{1}},
ImagePadding->{{40,1},{20,1}}]

(* equations of motion *)

DGL={t''[τ]==-a[t[τ]] ((x'[τ])^2) a'[t[τ]],
t'[0]==(Sqrt[1+(a0^2) ((ẋ0)^2)]),t[0]==t0 Gyr,
x''[τ]==-((2 t'[τ] x'[τ] a'[t[τ]])/(a[t[τ]])),
x'[0]==ẋ0, x[0]==0};

(* t[τ] comoving *)

ẋ0=0;sol1=NDSolve[DGL,{t,x},{τ,0,-t0 Gyr}];
Plot[(t[τ Gyr]/Gyr)/.sol1,{τ,0,-t0},
Frame->True,ImagePadding->{{40,1},{20,1}}]

(* t[τ] peculiar motion *)

ẋ0=999/1000;sol2=NDSolve[DGL,{t,x},{τ,0,-t0 Gyr}];
Plot[(t[τ Gyr]/Gyr)/.sol2,{τ,0,-2.0008223*^17/Gyr},
Frame->True,ImagePadding->{{40,1},{20,1}}]

(* x[τ] peculiar motion *)

ẋ0=999/1000;sol2=NDSolve[DGL,{t,x},{τ,0,-t0 Gyr}];
Plot[(-x[τ Gyr]a[t[τ Gyr]]/Pc)/.sol2,{τ,0,-2.0008223*^17/Gyr},
Frame->True,ImagePadding->{{40,1},{20,1}}]

Quit[]

(* flat DE+matter universe*)

ΩM=315/1000;
ΩΛ=1-ΩM;
H0=217618*^-23;
Gyr=(10^7)*36525*24*3600;
a0=1;

(* scale factor *)

a[t_]:=E^(H0 t);

(* t0 at a=1 *)

t0=-2.111*^17/Gyr;

(* scale factor plot *)

Plot[a[t Gyr],{t,-8t0,2t0},Frame->True,GridLines->{{t0},{1}}]

(* equations of motion *)

DGL={t''[τ]==-a[t[τ]] ((x'[τ])^2) a'[t[τ]],
t'[0]==(Sqrt[1+(a0^2) ((ẋ0)^2)]),t[0]==0,
x''[τ]==-((2 t'[τ] x'[τ] a'[t[τ]])/(a[t[τ]])),
x'[0]==ẋ0, x[0]==0};

(* t[τ] comoving *)

ẋ0=0;sol1=NDSolve[DGL,{t,x},{τ,0,8t0 Gyr}];
Plot[(t[τ Gyr]/Gyr)/.sol1,{τ,0,8t0},Frame->True]

(* t[τ] ultrarelativistiv *)

ẋ0=9999/10000;sol2=NDSolve[DGL,{t,x},{τ,0, -4*^17}];
Plot[(t[τ Gyr]/Gyr)/.sol2,{τ,0,-4*^17/Gyr },Frame->True]

ẋ0=99999/100000;sol2=NDSolve[DGL,{t,x},{τ,0,-4*^17}];
Plot[(x[τ Gyr]/Gyr)/.sol2,{τ,0,-4*^17/Gyr},Frame->True]
Plot[(x[τ Gyr] a[t[τ Gyr]])/.sol2,{τ,0,-4*^17/Gyr},Frame->True]

Plot[(a[t[τ Gyr]])/.sol2,{τ,0,-4*^17/Gyr},Frame->True]

Plot[(a[tt])/.sol2,{tt,0,-40*^17},Frame->True]