ord=10;
S=Table[0,{i,0,ord+2},{j,0,ord+2},{k,0,ord+2}];

Clear[i,j,k]

S[[3,1,1]]=-h[r]/2;
S[[1,3,1]]=f[r]/2;
S[[1,1,2]]=-p[r];

ST[i_,j_,k_]:=Sum[x*(i+2-x)*S[[x+1,y+1,z+1]]*S[[(i+2)-x+1,j-y+1,k-z+1]],{x,1,i+1},{y,0,j},{z,0,k}];
SR1[i_,j_,k_]:=Sum[D[S[[x+1,y+1,z+1]],r]*D[S[[i-x+1,j-y+1,k-z+1]],r],{x,0,i},{y,0,j},{z,0,k}];
SR2[i_,j_,k_]:=Sum[y*S[[x+1,y+1,z+1]]*D[S[[i-x+1,(j+1)-y+1,k-z+1]],r],{x,0,i},{y,1,j+1},{z,0,k}];
SR3[i_,j_,k_]:=Sum[y*(j+2-y)*S[[x+1,y+1,z+1]]*S[[i-x+1,(j+2)-y+1,k-z+1]],{x,0,i},{y,1,j+1},{z,0,k}];
SE1[i_,j_,k_]:=Sum[z*(k+1-z)*S[[x+1,y+1,z+1]]*S[[i-x+1,j-y+1,(k+1)-z+1]],{x,0,i},{y,0,j},{z,1,k}];
SE2[i_,j_,k_]:=Sum[z*(k-z)*S[[x+1,y+1,z+1]]*S[[i-x+1,j-y+1,k-z+1]],{x,0,i},{y,0,j},{z,1,k-1}];
div[i_,j_,k_]:=1-2*(i/2+j/2+k);
(*Coef[i_,j_,k_]:=If[(i>=2)||(j>=2)||(k>=1),(-1/(2*h[r])*ST[i,j,k]+1/(2*f[r])*(SR1[i,j,k]+2*SR2[i,j,k]+SR3[i,j,k])-1/(2*p[r])*(2*SE1[i,j,k]+SE2[i,j,k]))/div[i,j,k],S[[i+1,j+1,k+1]]];*)
Coef[i_,j_,k_]:=(-1/(2*h[r])*ST[i,j,k]+1/(2*f[r])*(SR1[i,j,k]+2*SR2[i,j,k]+SR3[i,j,k])-1/(2*p[r])*(2*SE1[i,j,k]+SE2[i,j,k]))/div[i,j,k];
(*S[[i+1,j+1,k+1]]=If[(i>2)||(j>2)||(k>1),(-1/(2*h[r])*ST+1/(2*f[r])*(SR1+2*SR2+SR3)-1/p[r]*(2*SE1+SE2))/div,S[[i+1,j+1,k+1]]]*)

ord=10;
S=Table[0,{i,0,ord+2},{j,0,ord+2},{k,0,ord+2}];
Clear[i,j,k];
ord=8;
S[[3,1,1]]=-h[r]/2;
S[[1,3,1]]=f[r]/2;
S[[1,1,2]]=-p[r];
For[l=3,l≤ord,l++,
For[i=0,i≤ord,i++,
For[j=0,j≤ord,j++,
For[k=0,k≤ord,k++,
S[[i+1,j+1,k+1]]=If[i+j+2*k==l,Coef[i,j,k],S[[i+1,j+1,k+1]]]//Simplify
];
];
];
Print[l];
];

S

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

SetAttributes[ϵ,Constant]

temp=-1/(2*h[r])*D[σ[t,r,η],t]^2+1/(2*f[r])*D[σ[t,r,η],r]^2-1/(2*p[r])*(D[σ[t,r,η],η]/ϵ^2)^2*(2(η*ϵ^2)+(η*ϵ^2)^2)

drule={D[dt[t/ϵ]*ϵ,t]->1,D[dr[r/ϵ]*ϵ,r]->1}

n=8;
Coefficient[σ[t,r,η],ϵ^n]/.drule//Simplify
Coefficient[temp,ϵ^n]/.drule//Simplify
Coefficient[temp-σ[t,r,η],ϵ^n]/.drule//Simplify

Save["Sigma.m", {S}]