(* 1. Set Parameters *)

$Assumptions=k>0;$Assumptions=η>0;
aE=1;HI=1;Γ=1/(((13824 (477087-97732 (Sqrt[21])))/8575)(HI^6));c2[k_]:=-I c1[k];

(* η=((3+(Sqrt[21]))/(6 aE HI))(τ^((9-(Sqrt[21]))/10)) *)

Plot[((3+(Sqrt[21]))/(6 aE HI))(τ^((9-(Sqrt[21]))/10)),{τ,0.6,1.2}]

(* 2. u0, u1 *)

FullSimplify[Expand[(solU1)-(u1[η,k] /. {c1[k]->c1,c2[k]->c2})]] (* check *)

u0[η_,k_]:=c1[k](Sqrt[k η])BesselJ[0,k η]+c2[k](Sqrt[k η])BesselY[0,k η];
u1[η_,k_]:=(((2^6)  ((((((Sqrt[21])-3)/2)))^9) π (HI^6))/(8575 (3^2)*7 ((((((Sqrt[21])-3)/2)aE HI η))^(9+(Sqrt[21])))))( c2 [k]((((k η))^(5/2)) BesselJ[0,k η] (-(213002167+46502521 (Sqrt[21])) HypergeometricPFQ[{(1/2),-(7/2)-((Sqrt[21])/2)},{1,1,-(5/2)-((Sqrt[21])/2)},-(k^2) (η^2)]-3 (34391077+7499743 (Sqrt[21])) HypergeometricPFQ[{(3/2),-(7/2)-((Sqrt[21])/2)},{2,2,-(5/2)-((Sqrt[21])/2)},-(k^2) (η^2)]+51450 (279+61 (Sqrt[21])) (k^2) (η^2) HypergeometricPFQ[{(1/2),-(5/2)-((Sqrt[21])/2)},{1,1,-(3/2)-((Sqrt[21])/2)},-(k^2) (η^2)]-85750 (819+179 (Sqrt[21])) (k^2) (η^2) HypergeometricPFQ[{(3/2),-(5/2)-((Sqrt[21])/2)},{2,2,-(3/2)-((Sqrt[21])/2)},-(k^2) (η^2)])+(1/(Sqrt[π]))(BesselJ[0,k η](-205800 (669+146 (Sqrt[21])) (k^4) (η^4)  MeijerG[{{(1/4) (15+2 (Sqrt[21]))},{(3/4),(3/4)}},{{(1/4),(1/4),(1/4)},{(3/4),(1/4) (11+2 (Sqrt[21]))}},k η,(1/2)]+2*171500 (3927+857 (Sqrt[21])) (k^3) (η^3)  MeijerG[{{(1/4) (17+2 (Sqrt[21]))},{(1/4),(1/4)}},{{-(1/4),-(1/4),(3/4)},{(1/4),(1/4) (13+2 (Sqrt[21]))}},k η,(1/2)]+14 (176254865+38465701 (Sqrt[21])) (k^2) (η^2)  MeijerG[{{(1/4) (19+2 (Sqrt[21]))},{(3/4),(3/4)}},{{(1/4),(1/4),(1/4)},{(3/4),(1/4) (15+2 (Sqrt[21]))}},k η,(1/2)]+84 (28445153+6206377 (Sqrt[21])) k η  MeijerG[{{(1/4) (21+2 (Sqrt[21]))},{(1/4),(1/4)}},{{-(1/4),-(1/4),(3/4)},{(1/4),(1/4) (17+2 (Sqrt[21]))}},k η,(1/2)])+BesselY[0,k η](102900 (669+146 (Sqrt[21])) (k^3) (η^3)  MeijerG[{{(5/4),(1/4) (17+2 (Sqrt[21]))},{(1/4)}},{{(3/4),(3/4)},{(1/4),(3/4),(1/4) (13+2 (Sqrt[21]))}},k η,(1/2)]-171500 (3927+857 (Sqrt[21])) (k^2) (η^2)  MeijerG[{{(3/4),(1/4) (19+2 (Sqrt[21]))},{-(1/4)}},{{(1/4),(5/4)},{-(1/4),(1/4),(1/4) (15+2 (Sqrt[21]))}},k η,(1/2)]-7 (176254865+38465701 (Sqrt[21])) k η MeijerG[{{(5/4),(1/4) (21+2 (Sqrt[21]))},{(1/4)}},{{(3/4),(3/4)},{(1/4),(3/4),(1/4) (17+2 (Sqrt[21]))}},k η,(1/2)]-42 (28445153+6206377 (Sqrt[21]))  MeijerG[{{(3/4),(1/4) (23+2 (Sqrt[21]))},{-(1/4)}},{{(1/4),(5/4)},{-(1/4),(1/4),(1/4) (19+2 (Sqrt[21]))}},k η,(1/2)])))+c1 [k]((((k η))^(5/2)) BesselY[0,k η] ((213002167+46502521 (Sqrt[21])) HypergeometricPFQ[{(1/2),-(7/2)-((Sqrt[21])/2)},{1,1,-(5/2)-((Sqrt[21])/2)},-(k^2) (η^2)]+3 (34391077+7499743 (Sqrt[21])) HypergeometricPFQ[{(3/2),-(7/2)-((Sqrt[21])/2)},{2,2,-(5/2)-((Sqrt[21])/2)},-(k^2) (η^2)]-51450 (279+61 (Sqrt[21])) (k^2) (η^2) HypergeometricPFQ[{(1/2),-(5/2)-((Sqrt[21])/2)},{1,1,-(3/2)-((Sqrt[21])/2)},-(k^2) (η^2)]+85750 (819+179 (Sqrt[21])) (k^2) (η^2) HypergeometricPFQ[{(3/2),-(5/2)-((Sqrt[21])/2)},{2,2,-(3/2)-((Sqrt[21])/2)},-(k^2) (η^2)])+(1/(Sqrt[π]))BesselJ[0,k η](-102900 (669+146 (Sqrt[21])) (k^3) (η^3)  MeijerG[{{(5/4),(1/4) (17+2 (Sqrt[21]))},{(1/4)}},{{(3/4),(3/4)},{(1/4),(3/4),(1/4) (13+2 (Sqrt[21]))}},k η,(1/2)]+171500 (3927+857 (Sqrt[21])) (k^2) (η^2)  MeijerG[{{(5/4),(1/4) (19+2 (Sqrt[21]))},{(3/4)}},{{(5/4),(5/4)},{(1/4),(1/4),(1/4) (15+2 (Sqrt[21]))}},k η,(1/2)]+7 (176254865+38465701 (Sqrt[21])) k η  MeijerG[{{(5/4),(1/4) (21+2 (Sqrt[21]))},{(1/4)}},{{(3/4),(3/4)},{(1/4),(3/4),(1/4) (17+2 (Sqrt[21]))}},k η,(1/2)]+42 (28445153+6206377 (Sqrt[21])) MeijerG[{{(5/4),(1/4) (23+2 (Sqrt[21]))},{(3/4)}},{{(5/4),(5/4)},{(1/4),(1/4),(1/4) (19+2 (Sqrt[21]))}},k η,(1/2)])));

(* 3. h0, h1 *)

h0[η_,k_]:=(aE^((3-(Sqrt[21]))/4))((((((Sqrt[21])-3)/2)aE HI η))^(-(1/2)))u0[η,k];
h1[η_,k_]:=(aE^((3-(Sqrt[21]))/4))((((((Sqrt[21])-3)/2)aE HI η))^(-(1/2)))u1[η,k];

(* 4. Power Spectrum *)

PS0[η_,k_]:=((((k/(aE HI))))^3)((Abs[h0[η,k]])^2);
PS[η_,k_]:=((((k/(aE HI))))^3)((Abs[h0[η,k]+Γ h1[η,k]])^2);

(* 5. Fig. 1 *)

Plot[Log[PS[ η,(E^(-1))]/PS[ 0.9,(E^(-1))]],{η,0.8,200}]

Plot[{Log[((PS[ η,(E^(-2))])/(PS[ 0.9,(E^(-2))]))],Log[((PS[ η,(E^(-0))])/(PS[ 0.9,(E^(-0))]))],Log[((PS[ η,(E^2)])/(PS[ 0.9,(E^2)]))]},{η,0.8,3}]

Plot[{Log[((PS[ η,(E^(-6))])/(PS[ 0.9,(E^(-6))]))],Log[((PS[ η,(E^(-4))])/(PS[ 0.9,(E^(-4))]))],Log[((PS[ η,(E^(-2))])/(PS[ 0.9,(E^(-2))]))]},{η,0.8,3}]

Fig1=Plot[{Log[PS[ η,(E^(-30))]/PS[ 0.9,(E^(-30))]],Log[PS[ η,(E^(-20))]/PS[ 0.9,(E^(-20))]],Log[PS[ η,(E^(-10))]/PS[ 0.9,(E^(-10))]]},{η,0.8,3}]

Export["TensorPt1.pdf",Fig1]

Export["TensorPt1.eps",Fig1]

(* 6. Fig.2 *)

Fig2=Plot[Log[PS[ 3.0,(E^n)]/PS[ 0.9,(E^n)]],{n,-30.1,-10}]

data=Table[{n,Log[PS[ 3.0,(E^n)]/PS[ 0.9,(E^n)]]},{n,-30,-10}]

ListPlot[data]

curve=Fit[data,{1,n,(n^2),(n^3)},n]

Plot[curve,{n,-30,-10}]

Export["TensorPt2.pdf",Fig2]

Export["TensorPt2.eps",Fig2]