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

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

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

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

<< 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;
];

U[t,r,η]

(*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/ϵ]*ϵ,t]->1,D[dr[r/ϵ]*ϵ,r]->1,D[dt[t/ϵ]*ϵ^2,{t,2}]->0,D[dr[r/ϵ]*ϵ^2,{r,2}]->0,D[dt[t/ϵ]*ϵ^3,{t,3}]->0,D[dr[r/ϵ]*ϵ^3,{r,3}]->0,D[dt[t/ϵ]*ϵ^4,{t,4}]->0,D[dr[r/ϵ]*ϵ^4,{r,4}]->0,D[dt[t/ϵ]*ϵ^5,{t,5}]->0,D[dr[r/ϵ]*ϵ^5,{r,5}]->0,D[dt[t/ϵ]*ϵ^6,{t,6}]->0,D[dr[r/ϵ]*ϵ^6,{r,6}]->0}
erule={ϵ->1}

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

xrule={D[x[t/ϵ,r/ϵ,η]*ϵ,t]->D[x2[t/ϵ,r/ϵ,η]*ϵ,t]/.drule//Expand,D[x[t/ϵ,r/ϵ,η]*ϵ^2,{t,2}]->D[x2[t/ϵ,r/ϵ,η]*ϵ^2,{t,2}]/.drule//Expand,D[x[t/ϵ,r/ϵ,η]*ϵ,r]->D[x2[t/ϵ,r/ϵ,η]*ϵ,r]/.drule//Expand,D[x[t/ϵ,r/ϵ,η]*ϵ^2,{r,2}]->D[x2[t/ϵ,r/ϵ,η]*ϵ^2,{r,2}]/.drule//Expand,D[x[t/ϵ,r/ϵ,η],η]->D[x2[t/ϵ,r/ϵ,η],η]/.drule//Expand,D[x[t/ϵ,r/ϵ,η],{η,2}]->D[x2[t/ϵ,r/ϵ,η],{η,2}]/.drule//Expand}
x2rule={x->x2}
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;*)

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

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

s[t,r,η]=x[t/ϵ,r/ϵ,η]*ϵ^2+σ3[t,r,η]//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/.b2rule/.brule//Simplify//Expand

n=-4;
Coefficient[KGG,ϵ,n]/.x2rule//Simplify//Expand