SetAttributes[ϵ,Constant];
SetAttributes[M,Constant];
SetAttributes[κ,Constant];
ord=2;

(*h[r_]:=1;
f[r_]:=1/b[r]^2;
p[r_]:=r^2/b[r];*)

b[r_]:=1/2*M/r

brule={M->(1-b)*r/2};

<< Sigma.m
S2=S;
Clear[S];

Clear[i,j,k]
For[i=0,i≤ord+4,i++,
For[j=0,j≤ord+4,j++,
For[k=0,k≤(ord+4)/2,k++,
Subscript[S,i,j,k]=S2[[i+1,j+1,k+1]];
];
];
];

σ[t,r,η]=Sum[If[i+j+2k≤ord+4,Subscript[S,i,j,k]*dt[t/ϵ]^i*dr[r/ϵ]^j*η^k*ϵ^(i+j+2*k),0],{i,0,ord+4},{j,0,ord+4},{k,0,(ord+4)/2}]//Expand;

σ2[t,r,η]=Series[σ[t,r,η],{ϵ,0,2}]//Normal//Expand

σ3[t,r,η]=σ[t,r,η]-σ2[t,r,η]//Expand;

<< Delta.m

Clear[i,j,k]
For[i=0,i≤ord+2,i++,
For[j=0,j≤ord+2,j++,
For[k=0,k≤(ord+2)/2,k++,
Subscript[c,i,j,k]=cs[[i+1,j+1,k+1]];
];
];
];

lDel[t,r,η]=Sum[If[i+j+2k≤ord+2,Subscript[c,i,j,k]*dt[t/ϵ]^i*dr[r/ϵ]^j*η^k*ϵ^(i+j+2*k),0],{i,0,ord+2},{j,0,ord+2},{k,0,(ord+2)/2}]//Expand;

Dim=Dimensions[lDel[t,r,η]][[1]]

uprule={dr[r/ϵ]->dr[r],dt[t/ϵ]->dt[t],x[t/ϵ,r/ϵ,η]->x[t,r,η]}
dnrule={dr[r]->dr[r/ϵ],dt[t]->dt[t/ϵ],x[t,r,η]->x[t/ϵ,r/ϵ,η]}

Clear[i];
Us=Table[0,{i,Dim}];
For[i=1,i≤Dim,i++,
Us[[i]]=Series[Exp[lDel[t,r,η][[i]]/.uprule],{ϵ,0,ord+2}]/.dnrule//Normal;
];

Clear[U]

U[t,r,η]=Us[[1]];
For[i=2,i≤Dim,i++,
temp=Series[Expand[U[t,r,η]*Us[[i]]],{ϵ,0,ord+2}]//Normal;
U[t,r,η]=temp//Expand;
];

(*Usig2[t,r,η]=Series[Expand[σ[t,r,η]^2*U[t,r,η]],{ϵ,0,ord+6}]//Normal//Expand;*)

prule={dr[r/ϵ]->dr[r/ϵ-r2/ϵ],dt[t/ϵ]->dt[t/ϵ-t2/ϵ]}
drule={D[dt[t/ϵ]*ϵ,t]->1,D[dr[r/ϵ]*ϵ,r]->1,D[dt[t/ϵ]*ϵ^2,{t,2}]->0,D[dr[r/ϵ]*ϵ^2,{r,2}]->0}
d2rule={D[dt[t/ϵ]*ϵ,t]->1,D[dr[r/ϵ]*ϵ,r]->1,D[dt[t/ϵ]*ϵ^2,{t,2}]->0,D[dr[r/ϵ]*ϵ^2,{r,2}]->0}
erule={ϵ->1}

x2[t/ϵ,r/ϵ,η]=σ2[t,r,η]/ϵ^2//Expand

xrule={D[x[t/ϵ,r/ϵ,η]*ϵ,t]->D[x2[t/ϵ,r/ϵ,η]*ϵ,t],D[x[t/ϵ,r/ϵ,η]*ϵ^2,{t,2}]->D[x2[t/ϵ,r/ϵ,η]*ϵ^2,{t,2}],D[x[t/ϵ,r/ϵ,η]*ϵ,r]->D[x2[t/ϵ,r/ϵ,η]*ϵ,r],D[x[t/ϵ,r/ϵ,η]*ϵ^2,{r,2}]->D[x2[t/ϵ,r/ϵ,η]*ϵ^2,{r,2}],D[x[t/ϵ,r/ϵ,η],η]->D[x2[t/ϵ,r/ϵ,η],η],D[x[t/ϵ,r/ϵ,η],{η,2}]->D[x2[t/ϵ,r/ϵ,η],{η,2}]}
xerule={x[t,r,η]->x}

R[r_]:=(4b[r]-4b[r]^2+4*r*b[r]*b'[r]-3*r^2*b'[r]^2+4*r^2*b[r]*b''[r])/(2*r^2)

R[r]//Expand

q[r_]:=h[r]*p[r]^2/f[r];

box[A_]:=-1/h[r]*D[A,{t,2}]+1/(2*h[r]*p[r]^2)*D[q[r],r]*D[A,r]+1/f[r]*D[A,{r,2}]-2*(η*ϵ^2+1)/p[r]*D[A,η]/ϵ^2-(2η*ϵ^2+(η*ϵ^2)^2)/p[r]*D[A,{η,2}]/ϵ^4;

dAdB[A_,B_]:=-1/h[r]*D[A,t]*D[B,t]+1/f[r]*D[A,r]*D[B,r]-(2*η*ϵ^2+η^2*ϵ^4)/p[r]*D[A,η]*D[B,η]/ϵ^4;

(*sig3G01[t,r,η]=Series[-1/(8*Pi^2)*R[t,r,η]/6*Usig2[t,r,η]/.drule//Expand,{ϵ,0,ord+4}]//Normal//Expand;*)

(*sig3G02[t,r,η]=1/(8*Pi^2)*Series[box[Usig2[t,r,η]]/.drule//Expand,{ϵ,0,ord+4}]//Normal//Expand;*)

(*sig3G03[t,r,η]=Series[-3/(8*Pi^2)*U[t,r,η]*σ[t,r,η]*box[σ[t,r,η]]/.drule//Expand,{ϵ,0,ord+4}]//Normal//Expand;*)

(*sig3G04[t,r,η]=Series[-3/(4*Pi^2)*dAdB[U[t,r,η],σ[t,r,η]]*σ[t,r,η]/.drule//Expand,{ϵ,0,ord+4}]//Normal//Expand;*)

(*sig3G0[t,r,η]=sig3G01[t,r,η]+sig3G02[t,r,η]+sig3G03[t,r,η]+sig3G04[t,r,η];*)

b[r_]:=1-2*M/r;

(*Clear[b]*)

s[t,r,η]=x[t/ϵ,r/ϵ,η]*ϵ^2+σ3[t,r,η]//Expand;

(*G0[t,r,η]=Series[1/(8*Pi^2)*U[t,r,η]/s[t,r,η]/.uprule//Expand,{ϵ,0,ord}]/.dnrule//Normal//Expand;

(*n=2;
Coefficient[G0[t,r,η],ϵ,n]//Expand

(*G2[t,r,η]=Series[1/(8*Pi^2)*κ^4/180*(box[U[t,r,η]*(2ϵ^2*dt[r/ϵ]^2-σ[t,r,η])]-R[t,r,η]/6*U[t,r,η]*(2ϵ^2*dt[r/ϵ]^2-σ[t,r,η]))/.drule//Expand,{ϵ,0,ord}]//Normal//Expand;

Grn[t,r,η]=Series[1/(8*Pi^2)*U[t,r,η]*((*1/s[t,r,η]+*)κ^2/6(*+κ^4/180*(-2*ϵ^2*dt[t/ϵ]^2-s[t,r,η])*))/.brule/.uprule//Expand,{ϵ,0,2}]/.dnrule//Normal//Expand

KGG=box[Grn[t,r,η]]-R[r]*Grn[t,r,η]/6/.xrule/.drule//Expand

n=1;
Coefficient[KGG,ϵ,n]/.brule/.erule/.xerule//Expand//Simplify