SetAttributes[ϵ,Constant];
SetAttributes[M,Constant];
SetAttributes[κ,Constant];
ordn=-4;
ordk=2;

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

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

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

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

σ[t,r,η]=Sum[If[i+j+2k≤ordn+10,Subscript[S,i,j,k]*dt[t/ϵ]^i*dr[r/ϵ]^j*η^k*ϵ^(i+j+2*k),0],{i,0,ordn+10},{j,0,ordn+10},{k,0,(ordn+10)/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≤ordn+8,i++,
For[j=0,j≤ordn+8,j++,
For[k=0,k≤(ordn+8)/2,k++,
Subscript[c,i,j,k]=cs[[i+1,j+1,k+1]];
];
];
];

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

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

uprule={dr[r/ϵ]->dr[r],dt[t/ϵ]->dt[t],dr[r/ϵ-r2/ϵ]->dr[r-r2],dt[t/ϵ-t2/ϵ]->dt[t-t2],x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]->x[t,t2,r,r2,η],A_[(θ/ϵ),((θ2)/ϵ),(ϕ/ϵ),((ϕ2)/ϵ)]->A[θ,θ2,ϕ,ϕ2]}
dnrule={dr[r]->dr[r/ϵ],dt[t]->dt[t/ϵ],dr[r-r2]->dr[r/ϵ-r2/ϵ],dt[t-t2]->dt[t/ϵ-t2/ϵ],x[t,t2,r,r2,η]->x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η],A_[θ,θ2,ϕ,ϕ2]->A[(θ/ϵ),((θ2)/ϵ),(ϕ/ϵ),((ϕ2)/ϵ)]}

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

Clear[U]

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

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

prule={dr[r/ϵ]->dr[r/ϵ-r2/ϵ],dt[t/ϵ]->dt[t/ϵ-t2/ϵ]}
(*drule={D[dt[t/ϵ-t2/ϵ]*ϵ,t]->1,D[dr[r/ϵ-r2/ϵ]*ϵ,r]->1,D[dt[t/ϵ-t2/ϵ]*ϵ^2,{t,2}]->0,D[dr[r/ϵ-r2/ϵ]*ϵ^2,{r,2}]->0,D[dt[t/ϵ-t2/ϵ]*ϵ,t2]->-1,D[dr[r/ϵ-r2/ϵ]*ϵ,r2]->-1,D[dt[t/ϵ-t2/ϵ]*ϵ^2,{t2,2}]->0,D[dr[r/ϵ-r2/ϵ]*ϵ^2,{r2,2}]->0}*)
drule={D[dt[t/ϵ-t2/ϵ]*ϵ,t]->1,D[dr[r/ϵ-r2/ϵ]*ϵ,r]->1,D[dt[t/ϵ-t2/ϵ]*ϵ^2,{t,2}]->0,D[dr[r/ϵ-r2/ϵ]*ϵ^2,{r,2}]->0,D[dt[t/ϵ-t2/ϵ]*ϵ^3,{t,3}]->0,D[dr[r/ϵ-r2/ϵ]*ϵ^3,{r,3}]->0,D[dt[t/ϵ-t2/ϵ]*ϵ^4,{t,4}]->0,D[dr[r/ϵ-r2/ϵ]*ϵ^4,{r,4}]->0,
D[dt[t/ϵ-t2/ϵ]*ϵ^5,{t,5}]->0,D[dr[r/ϵ-r2/ϵ]*ϵ^5,{r,5}]->0,
D[dt[t/ϵ-t2/ϵ]*ϵ^6,{t,6}]->0,D[dr[r/ϵ-r2/ϵ]*ϵ^6,{r,6}]->0}
erule={ϵ->1}
r2rule={r2->r-ϵ*dr[r/ϵ]}
aprule={dr[r/ϵ-r2/ϵ]->dr[r/ϵ],dt[t/ϵ-t2/ϵ]->dt[t/ϵ],dt[(t-t2)/ϵ]->dt[t/ϵ],x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]->x[t/ϵ,r/ϵ,η]}
up2rule={dr[r/ϵ]->dr[r],dt[t/ϵ]->dt[t],x[t/ϵ,r/ϵ,η]->x[t,r,η],A_[(θ/ϵ),((θ2)/ϵ),(ϕ/ϵ),((ϕ2)/ϵ)]->A[θ,θ2,ϕ,ϕ2]}
dn2rule={dr[r]->dr[r/ϵ],dt[t]->dt[t/ϵ],x[t,r,η]->x[t/ϵ,r/ϵ,η],A_[θ,θ2,ϕ,ϕ2]->A[(θ/ϵ),((θ2)/ϵ),(ϕ/ϵ),((ϕ2)/ϵ)]}

x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]=σ2[t,r,η]/ϵ^2/.prule;

xrule={D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,t]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,t],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{t,2}]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{t,2}],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,t2]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,t2],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{t2,2}]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{t2,2}],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,r]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,r],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{r,2}]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{r,2}],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,r2]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ,r2],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{r2,2}]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2,{r2,2}],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η],η]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η],η],D[x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η],{η,2}]:>D[x2[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η],{η,2}]}

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

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;

box2[A_]:=-D[A,{t2,2}]+b[r2]^2*D[A,{r2,2}]+2*b[r2]^2/r2*D[A,r2]-2*(η*ϵ^2+1)*b[r2]/r2^2*D[A,η]/ϵ^2-(2η*ϵ^2+(η*ϵ^2)^2)*b[r2]/r2^2*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;

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

boxf[A_]:=-D[A,{t,2}]+b[r]^2*D[A,{r,2}]+2*b[r]^2/r*D[A,r]+b[r]/r^2*D[A/.y1rule,{θ,2}]+b[r]/r^2*Cot[θ]*D[A/.y1rule,θ]+b[r]/(r^2*Sin[θ]^2)*D[A/.y1rule,{ϕ,2}]/.y2rule;

box2f[A_]:=-D[A,{t2,2}]+b[r2]^2*D[A,{r2,2}]+2*b[r2]^2/r2*D[A,r2]+b[r2]/r2^2*D[A/.y1rule,{θ2,2}]+b[r2]/r2^2*Cot[θ2]*D[A/.y1rule,θ2]+b[r2]/(r2^2*Sin[θ2]^2)*D[A/.y1rule,{ϕ2,2}]/.y2rule;

dAdBf[A_,B_]:=-D[A,t]*D[B,t]+b[r]^2*D[A,r]*D[B,r]+b[r]/r^2*D[A/.y1rule,θ]*D[B/.y1rule,θ]+b[r]/(r^2*Sin[θ]^2)*D[A/.y1rule,ϕ]*D[B/.y1rule,ϕ]/.y2rule;

dAdB2f[A_,B_]:=-D[A,t2]*D[B,t2]+b[r2]^2*D[A,r2]*D[B,r2]+b[r2]/r2^2*D[A/.y1rule,θ2]*D[B/.y1rule,θ2]+b[r2]/(r2^2*Sin[θ2]^2)*D[A/.y1rule,ϕ2]*D[B/.y1rule,ϕ2]/.y2rule;

Mp[A_,B_,s_,k_]:=Sum[(Subscript[A, i,j])*(Subscript[B, s-i,k-j]),{i,-8,ordn+6},{j,0,k}]//Expand;

u[t_,t2_,r_,r2_,η_]:=Series[U[t,r,η]/.prule,{ϵ,0,ordn+8}]//Normal//Expand;

s[t_,t2_,r_,r2_,η_]:=x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2+σ3[t,r,η]/.prule;

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

y1rule={η->Y[θ/ϵ,θ2/ϵ,ϕ/ϵ,ϕ2/ϵ]};
y2rule={Y[θ/ϵ,θ2/ϵ,ϕ/ϵ,ϕ2/ϵ]->η};

bGrn=Series[box[Grn]/.xrule/.drule/.uprule//Expand,{ϵ,0,ordn+4}]/.dnrule//Normal//Expand;
bpGrn=Series[box2[Grn]/.xrule/.drule/.uprule//Expand,{ϵ,0,ordn+4}]/.dnrule//Normal//Expand;

bbpGrn=Series[box[bpGrn]/.xrule/.drule/.uprule//Expand,{ϵ,0,ordn+2}]/.dnrule//Normal//Expand;

Clear[i,j,k];
For[i=1,i≤4,i++,
For[j=1,j≤4,j++,
For[k=1,k≤4,k++,
(Subscript[A, i,j,k])=0;
(Subscript[Ap, i,j,k])=0;
];
(Subscript[R, i,j])=0;
(Subscript[Rp, i,j])=0;
(Subscript[g, i,j])=0;
(Subscript[gp, i,j])=0;
];
];
(Subscript[A, 1,3,3])=3*M*Sin[θ]^2 - r*Sin[θ]^2;
(Subscript[A, 3,2,3])=Cot[θ];
(Subscript[A, 3,3,2])=Cot[θ];
(Subscript[A, 1,2,2])=3*M - r;
(Subscript[A, 3,3,1])=(3*M - r)/((2*M - r)*r);
(Subscript[A, 3,1,3])=(3*M - r)/((2*M - r)*r);
(Subscript[A, 2,3,3])=-(Cos[θ]*Sin[θ]);
(Subscript[A, 2,2,1])=(3*M - r)/((2*M - r)*r);
(Subscript[A, 2,1,2])=(3*M - r)/((2*M - r)*r);
(Subscript[A, 1,1,1])=(2*M)/((2*M - r)*r);

(Subscript[Ap, 1,3,3])=3*M*Sin[θ2]^2 - r2*Sin[θ2]^2;
(Subscript[Ap, 3,2,3])=Cot[θ2];
(Subscript[Ap, 3,3,2])=Cot[θ2];
(Subscript[Ap, 1,2,2])=3*M - r2;
(Subscript[Ap, 3,3,1])=(3*M - r2)/((2*M - r2)*r2);
(Subscript[Ap, 3,1,3])=(3*M - r2)/((2*M - r2)*r2);
(Subscript[Ap, 2,3,3])=-(Cos[θ2]*Sin[θ2]);
(Subscript[Ap, 2,2,1])=(3*M - r2)/((2*M - r2)*r2);
(Subscript[Ap, 2,1,2])=(3*M - r2)/((2*M - r2)*r2);
(Subscript[Ap, 1,1,1])=(2*M)/((2*M - r2)*r2);

(Subscript[R, 1,1])=(2*(3*M^2 - 2*M*r))/((2*M - r)^2*r^2);
(Subscript[R, 2,2])=(2*(3*M^2 - M*r))/((2*M - r)*r);
(Subscript[R, 3,3])=(2*(3*M^2*Sin[θ]^2 - M*r*Sin[θ]^2))/((2*M - r)*r);

(Subscript[Rp, 1,1])=(2*(3*M^2 - 2*M*r2))/((2*M - r2)^2*r2^2);
(Subscript[Rp, 2,2])=(2*(3*M^2 - M*r2))/((2*M - r2)*r2);
(Subscript[Rp, 3,3])=(2*(3*M^2*Sin[θ2]^2 - M*r2*Sin[θ2]^2))/((2*M - r2)*r2);

(Subscript[g, 1,1])=1/b[r]^2
(Subscript[g, 2,2])=r^2/b[r];
(Subscript[g, 3,3])=r^2*Sin[θ]^2/b[r];
(Subscript[g, 4,4])=-1;

(Subscript[gp, 1,1])=1/b[r2]^2
(Subscript[gp, 2,2])=r2^2/b[r2];
(Subscript[gp, 3,3])=r2^2*Sin[θ2]^2/b[r2];
(Subscript[gp, 4,4])=-1;

N4=Table[0,{i,10}];
N2=Table[0,{i,11},{j,4},{k,4}];
N2p=Table[0,{i,11},{j,4},{k,4}];
N0=Table[0,{i,13},{j,4},{k,4},{j2,4},{k2,4}];
NKcomp=Table[0,{j,4},{k,4},{j2,4},{k2,4}];

clist={{4,4,4,4},{1,1,4,4},{2,2,4,4},{3,3,4,4}}
cnum=4

c=clist[[cnum]]

c1={r,θ,ϕ,t};
c2={r2,θ2,ϕ2,t2};

(Subscript[v, 1]):=c1[[c[[1]]]];
(Subscript[v, 2]):=c1[[c[[2]]]];
(Subscript[v, 3]):=c2[[c[[3]]]];
(Subscript[v, 4]):=c2[[c[[4]]]];

(*Initialization*)
For[i=-20,i≤20,i++,
For[j=0,j≤4,j++,
(Subscript[G, i,j])=0;
(Subscript[G1000, i,j])=0;
(Subscript[G0100, i,j])=0;
(Subscript[G1100, i,j])=0;
(Subscript[G0010, i,j])=0;
(Subscript[G0001, i,j])=0;
(Subscript[G0011, i,j])=0;
(Subscript[G1010, i,j])=0;
(Subscript[G1001, i,j])=0;
(Subscript[G1011, i,j])=0;
(Subscript[G0110, i,j])=0;
(Subscript[G0101, i,j])=0;
(Subscript[G0111, i,j])=0;
(Subscript[G1110, i,j])=0;
(Subscript[G1101, i,j])=0;
(Subscript[G1111, i,j])=0;
(Subscript[Gsq, i,j])=0;
(Subscript[bG, i,j])=0;
(Subscript[bpG, i,j])=0;
(Subscript[bbpG, i,j])=0;
];
];

For[i=-2,i≤ordn+6,i++,
For[j=0,j≤4,j++,
(Subscript[G, i,j])=Coefficient[Coefficient[Grn,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(*Print[(Subscript[G, i,j])];*)
];
];

(*Note: Need to deal with r2 order below!*)

Grn1000=D[Grn/.y1rule,(Subscript[v, 1])]/.y2rule/.xrule/.drule//Expand;
Grn0100=D[Grn/.y1rule,(Subscript[v, 2])]/.y2rule/.xrule/.drule//Expand;
Grn1100=(D[Grn1000/.y1rule,(Subscript[v, 2])]-Sum[(Subscript[A, i,c[[1]],c[[2]]])*D[Grn/.y1rule,c1[[i]]],{i,1,4}])/.y2rule/.xrule/.drule//Expand;
Grn0010=D[Grn/.y1rule,(Subscript[v, 3])]/.y2rule/.xrule/.drule//Expand;
Grn0001=D[Grn/.y1rule,(Subscript[v, 4])]/.y2rule/.xrule/.drule//Expand;
Grn0011=(D[Grn0010/.y1rule,(Subscript[v, 4])]-Sum[(Subscript[Ap, i,c[[3]],c[[4]]])*D[Grn/.y1rule,c2[[i]]],{i,1,4}])/.y2rule/.xrule/.drule//Expand;
Grn1010=D[Grn1000/.y1rule,(Subscript[v, 3])]/.y2rule/.xrule/.drule//Expand;
Grn1001=D[Grn1000/.y1rule,(Subscript[v, 4])]/.y2rule/.xrule/.drule//Expand;
Grn1011=(D[Grn1010/.y1rule,(Subscript[v, 4])]-Sum[(Subscript[Ap, i,c[[1]],c[[2]]])*D[Grn1000/.y1rule,c2[[i]]],{i,1,4}])/.y2rule/.xrule/.drule//Expand;
Grn0110=D[Grn0100/.y1rule,(Subscript[v, 3])]/.y2rule/.xrule/.drule//Expand;
Grn0101=D[Grn0100/.y1rule,(Subscript[v, 4])]/.y2rule/.xrule/.drule//Expand;
Grn0111=(D[Grn0110/.y1rule,(Subscript[v, 4])]-Sum[(Subscript[Ap, i,c[[1]],c[[2]]])*D[Grn0100/.y1rule,c2[[i]]],{i,1,4}])/.y2rule/.xrule/.drule//Expand;
Grn1110=D[Grn1100/.y1rule,(Subscript[v, 3])]/.y2rule/.xrule/.drule//Expand;
Grn1101=D[Grn1100/.y1rule,(Subscript[v, 4])]/.y2rule/.xrule/.drule//Expand;
Grn1111=(D[Grn1110/.y1rule,(Subscript[v, 4])]-Sum[(Subscript[Ap, i,c[[1]],c[[2]]])*D[Grn1100/.y1rule,c2[[i]]],{i,1,4}])/.y2rule/.xrule/.drule//Expand;

(*Derivatives*)
For[i=-8,i≤ordn+6,i++,
For[j=0,j≤4,j++,
(Subscript[G1000, i,j])=Coefficient[Coefficient[Grn1000,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G0100, i,j])=Coefficient[Coefficient[Grn0100,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G1100, i,j])=Coefficient[Coefficient[Grn1100,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G0010, i,j])=Coefficient[Coefficient[Grn0010,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G0001, i,j])=Coefficient[Coefficient[Grn0001,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G0011, i,j])=Coefficient[Coefficient[Grn0011,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G1010, i,j])=Coefficient[Coefficient[Grn1010,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G1001, i,j])=Coefficient[Coefficient[Grn1001,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G1011, i,j])=Coefficient[Coefficient[Grn1011,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G0110, i,j])=Coefficient[Coefficient[Grn0110,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G0101, i,j])=Coefficient[Coefficient[Grn0101,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G0111, i,j])=Coefficient[Coefficient[Grn0111,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G1110, i,j])=Coefficient[Coefficient[Grn1110,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G1101, i,j])=Coefficient[Coefficient[Grn1101,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[G1111, i,j])=Coefficient[Coefficient[Grn1111,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
];
];

(Subscript[G0011, -3,0])

D[1/(x[t/ϵ,t2/ϵ,r/ϵ,r2/ϵ,η]*ϵ^2),t]/.xrule/.drule

For[i=-4,i≤ordn+4,i++,
For[j=0,j≤4,j++,
(Subscript[Gsq, i,j])=Mp[G,G,i,j];
];
];

For[i=-8,i≤ordn+4,i++,
For[j=0,j≤4,j++,
(Subscript[bG, i,j])=Coefficient[Coefficient[bGrn,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[bpG, i,j])=Coefficient[Coefficient[bpGrn,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
(Subscript[bbpG, i,j])=Coefficient[Coefficient[bbpGrn,κ,j],ϵ,i]*κ^j*ϵ^i//Expand;
];
];

GsqS=Sum[(Subscript[Gsq, i,ordk]),{i,-8,0}];
N4[[1]]=Series[1/2*box[box2[GsqS]/.xrule/.drule]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[2]]=Series[3*Sum[(Subscript[bG, i,j])*(Subscript[bpG, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[3]]=Series[3*Sum[(Subscript[G, i,j])*(Subscript[bbpG, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[4]]=Series[-6*Sum[dAdB2[(Subscript[bG, i,j]),(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[5]]=Series[-6*Sum[dAdB[(Subscript[bpG, i,j]),(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[6]]=Series[R[r]*Sum[dAdB2[(Subscript[G, i,j]),(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[7]]=Series[R[r2]*Sum[dAdB[(Subscript[G, i,j]),(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[8]]=Series[-2*R[r]*Sum[(Subscript[G, i,j])*(Subscript[bpG, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[9]]=Series[-2*R[r2]*Sum[(Subscript[G, i,j])*(Subscript[bG, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N4[[10]]=Series[1/2*R[r]*R[r2]*GsqS/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[1,c[[1]],c[[2]]]]=Series[-8*Sum[dAdB2f[(Subscript[G0100, i,j]),(Subscript[G1000, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[1,c[[3]],c[[4]]]]=Series[-8*Sum[dAdBf[(Subscript[G0010, i,j]),(Subscript[G0001, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[2,c[[1]],c[[2]]]]=Series[8*Sum[(Subscript[G0100, i,j])*box2f[(Subscript[G1000, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[2,c[[3]],c[[4]]]]=Series[8*Sum[(Subscript[G0001, i,j])*boxf[(Subscript[G0010, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[3,c[[1]],c[[2]]]]=Series[8*Sum[(Subscript[G1000, i,j])*box2f[(Subscript[G0100, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[3,c[[3]],c[[4]]]]=Series[8*Sum[(Subscript[G0010, i,j])*boxf[(Subscript[G0001, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[4,c[[1]],c[[2]]]]=Series[4*Sum[dAdB2f[(Subscript[G, i,j]),(Subscript[G1100, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[4,c[[3]],c[[4]]]]=Series[4*Sum[dAdBf[(Subscript[G, i,j]),(Subscript[G0011, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

(*N2[[5,c[[1]],c[[2]]]]=Series[-4*Sum[box2[(Subscript[G, i,j])]*(Subscript[G1100, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.uprule//Expand,{ϵ,0,ordn}]/.dnrule//Normal//Expand*)
N2[[5,c[[1]],c[[2]]]]=Series[-4*Sum[(Subscript[bpG, i,j])*(Subscript[G1100, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

(*N2p[[5,c[[3]],c[[4]]]]=Series[-4*Sum[box[(Subscript[G, i,j])]*(Subscript[G0011, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.uprule//Expand,{ϵ,0,ordn}]/.dnrule//Normal//Expand*)
N2p[[5,c[[3]],c[[4]]]]=Series[-4*Sum[(Subscript[bG, i,j])*(Subscript[G0011, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[6,c[[1]],c[[2]]]]=Series[-4*Sum[(Subscript[G, i,j])*box2f[(Subscript[G1100, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[6,c[[3]],c[[4]]]]=Series[-4*Sum[(Subscript[G, i,j])*boxf[(Subscript[G0011, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[7,c[[1]],c[[2]]]]=Series[-4*R[r2]*Sum[(Subscript[G1000, i,j])*(Subscript[G0100, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[7,c[[3]],c[[4]]]]=Series[-4*R[r]*Sum[(Subscript[G0010, i,j])*(Subscript[G0001, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[8,c[[1]],c[[2]]]]=Series[2*R[r2]*Sum[(Subscript[G, i,j])*(Subscript[G1100, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[8,c[[3]],c[[4]]]]=Series[2*R[r]*Sum[(Subscript[G, i,j])*(Subscript[G0011, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[9,c[[1]],c[[2]]]]=Series[-2*(Subscript[R, c[[1]],c[[2]]])*Sum[dAdB2[(Subscript[G, i,j]),(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[9,c[[3]],c[[4]]]]=Series[-2*(Subscript[Rp, c[[3]],c[[4]]])*Sum[dAdB[(Subscript[G, i,j]),(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[10,c[[1]],c[[2]]]]=Series[4*(Subscript[R, c[[1]],c[[2]]])*Sum[(Subscript[G, i,j])*box2[(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[10,c[[3]],c[[4]]]]=Series[4*(Subscript[Rp, c[[3]],c[[4]]])*Sum[(Subscript[G, i,j])*box[(Subscript[G, k-i,ordk-j])],{i,-8,ordn+6},{j,0,ordk},{k,-6,-2}]
/.xrule/.drule/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2[[11,c[[1]],c[[2]]]]=Series[-R[r2]*(Subscript[R, c[[1]],c[[2]]])*GsqS/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N2p[[11,c[[3]],c[[4]]]]=Series[-R[r]*(Subscript[Rp, c[[3]],c[[4]]])*GsqS/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[1,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[4*Sum[(Subscript[G0110, i,j])*(Subscript[G1001, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[2,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[4*Sum[(Subscript[G1010, i,j])*(Subscript[G0101, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[3,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[Sum[(Subscript[G1100, i,j])*(Subscript[G0011, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[4,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[Sum[(Subscript[G, i,j])*(Subscript[G1111, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[5,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[Sum[-2*(Subscript[G0100, i,j])*(Subscript[G1011, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[6,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[Sum[-2*(Subscript[G1000, i,j])*(Subscript[G0111, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[7,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[Sum[-2*(Subscript[G0001, i,j])*(Subscript[G1110, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[8,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[Sum[-2*(Subscript[G0010, i,j])*(Subscript[G1101, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[9,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[2*(Subscript[Rp, c[[3]],c[[4]]])*Sum[(Subscript[G1000, i,j])*(Subscript[G0100, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[10,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[2*(Subscript[R, c[[1]],c[[2]]])*Sum[(Subscript[G0010, i,j])*(Subscript[G0001, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[11,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[-(Subscript[Rp, c[[3]],c[[4]]])*Sum[(Subscript[G1100, i,j])*(Subscript[G, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[12,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[-(Subscript[R, c[[1]],c[[2]]])*Sum[(Subscript[G0011, i,j])*(Subscript[G, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

N0[[13,c[[1]],c[[2]],c[[3]],c[[4]]]]=Series[(Subscript[Rp, c[[3]],c[[4]]])*(Subscript[R, c[[1]],c[[2]]])*Sum[(Subscript[G, i,j])*(Subscript[G, k-i,ordk-j]),{i,-8,ordn+6},{j,0,ordk},{k,-8,-4}]/.aprule/.r2rule/.up2rule,{ϵ,0,ordn}]/.dn2rule//Normal//Expand

Clear[i]
NKcomp[[c[[1]],c[[2]],c[[3]],c[[4]]]]=Sum[N0[[i,c[[1]],c[[2]],c[[3]],c[[4]]]]/72,{i,1,13}]+(Subscript[gp, c[[3]],c[[4]]])*Sum[N2[[i,c[[1]],c[[2]]]]/288,{i,1,11}]+(Subscript[g, c[[1]],c[[2]]])*Sum[N2p[[i,c[[3]],c[[4]]]]/288,{i,1,11}]+(Subscript[g, c[[1]],c[[2]]])*(Subscript[gp, c[[3]],c[[4]]])*Sum[N4[[i]]/288,{i,1,10}]/.aprule/.r2rule/.up2rule/.uprule//Expand

(*x2rule={x[t,r,η]->σ2[t,r,η]/.erule}*)

Save["NKcompK"<>ToString[ordk]<>".m",{N4, N2, N2p, N0,NKcomp}];