Simplifying the Einstein tensor in case of a perturbed FRW metric

Question

I use the code in this thread's answer:

(Calculating Einstein tensor components in Kaluza-Klein model) to get the Einstein tensor components of a four-dimensional Kaluza Klein model. But instead of the diagonal matrix in that thread, I use a perturbed metric 0209156:

$ g_{00}= a^2(t) ( 1+ 2 \phi^{(1)} + \phi^{(2)} )$

$ g_{0i}=0$

$g_{ij}= a^2(t) [( 1- 2 \psi^{(1)}- \psi^{(2)})~\delta_{ij} + \frac{1}{2} D_{ij} \chi_{ij}^{(2)}] $,

Where:

$ D_{ij} =\partial_i \partial_j - \frac{1}{3} \delta_{ij} \partial^k \partial_k $

Here is the entire code I use:

coordList = {t, x, w, z};
coord[i_] := coordList[[i]]
met = a[t]^2 Exp[y] {{-(1 + 2 n [CapitalPhi][t, x, w, z, y]), 0, 0, 0},
    {0, 1 - 2 n [CapitalPsi][t, x, w, z, y] +  n 1/2 D[Xi[t, x, w, z, y], {x, 2}] - 
      1/2 n D[Xi[t, x, w, z, y], {w, 2}] - 1/2 n D[Xi[t, x, w, z, y], {z, 2}], 
     n D[D[Xi[t, x, w, z, y], x], w], n D[D[Xi[t, x, w, z, y], x], z]}, {0, n D[D[Xi[t, x, w, z, y], x], w], 1 - 2 n [CapitalPsi][t, x, w, z, y] + n 1/2 D[Xi[t, x, w, z, y], {x, 2}] -  1/2 n D[Xi[t, x, w, z, y], {w, 2}] -  1/2 n D[Xi[t, x, w, z, y], {z, 2}], 
     n D[D[Xi[t, x, w, z, y], w], z]},{0, n D[D[Xi[t, x, w, z, y], x], z],n D[D[Xi[t, x, w, z, y], w], z],  1 - 2 n [CapitalPsi][t, x, w, z, y] + n 1/2 D[X[t, x, w, z, y], {x, 2}] - 
      1/2 n D[Xi[t, x, w, z, y], {w, 2}] -1/2 n D[Xi[t, x, w, z, y], {z, 2}]}};
invmetric = Inverse@metric;
m[a_, b_] := metric[[a, b]];
im[a_, b_] := invmetric[[a, b]];
d[a_, f_] := D[f, coord[a]];
dDel[b_, f_] := covd[f, coord[b]];
covd[0, x_] := 0;
d4[f_] := D[f, y]
term1[a_, b_] := dDel[b, d[a, f[t, y]]]/f[t, y];
term2[a_, b_] := d4[f[t, y]] d4[m[a, b]]/f[t, y];
term3[a_, b_] := -d4[d4[m[a, b]]];
term4[a_, b_, g_, d_] := im[g, d] (d4[m[a, g] d4[m[b, d]]]);
term5[a_, b_, g_, d_] := -im[g, d] (d4[m[g, d] d4[m[a, b]]])/2;
term6[a_, b_, g_, d_] := 
  m[a, b] (d4[im[g, d]] d4[m[g, d]] + (im[g, d] d4[m[g, d]])^2)/4;
ein[a_, b_] := -(term1[a, b] + term2[a, b] + term3[a, b] + 
       Sum[term4[a, b, g, d] + term5[a, b, g, d] + 

Where n is the order of perturbation. The problem now the output is huge. How to linearize such output by Series[Exp[n], {n,0,1}], and how to divide it to first order and second order perturbation terms like for instance in 0209156 Equations (30-40)

Edit

Here is an example of a code calculating Einstein tensor from the christoffel symbols, where it uses Series[--, {n,0,1}] for linearizton. And I do not get if it can divide the first order and second order perturbation terms. The idea I can not use this code because I substitute directly by the metric to the Einstein tensor equation. So this code needs to be modified for my calculation.

    x = {t, r, \[Theta], \[Phi]};
    met ={};
    metI=Inverse[met];
 Chrest = 
  ParallelTable[
   Normal[Series[
     1/2 (Sum[
        metI[[i, 
          p]] (D[met[[p, j]], x[[k]]] + D[met[[k, p]], x[[j]]] - 
           D[met[[j, k]], x[[p]]]), {p, 1, 4}]), {n, 0, 1}]], {i,1,4}, {j, 1, 4}, {k, 1, 4}];

Reim = ParallelTable[
   Normal[Series[
     D[Chrest[[i, l, j]], x[[k]]] - D[Chrest[[i, k, j]], x[[l]]] + 
      Sum[Chrest[[i, k, m]] Chrest[[m, l, j]], {m, 1, 4}] - 
      Sum[Chrest[[i, l, m]] Chrest[[m, k, j]], {m, 1, 4}], {n, 0, 
      1}]], {i, 1, 4}, {j, 1, 4}, {k, 1, 4}, {l, 1, 4}];

 Ricci =  ParallelTable[
   Normal[Series[Sum[Reim[[m, i, m, j]], {m, 1, 4}], {n, 0, 1}]], {i,1, 4}, {j, 1, 4}];

RScal =  Normal[Series[
    Sum[Sum[metI[[i, j]] Ricci[[i, j]], {i, 1, 4}], {j, 1, 4}], {n,0,1}]];

 Eins =   ParallelTable[
   Normal[Series[
     Ricci[[i, j]] - 1/2 met[[i, j]] RScal, {n, 0, 1}]], {i, 1, 
    4}, {j, 1, 4}];


EinsUD = ParallelTable[
       Normal[Series[
         Sum[metI[[p, i]] Eins[[p, j]], {p, 1, 4}], {n, 0, 1}]], {i, 1, 
        4}, {j, 1, 4}];
For[i = 1, i < 5, i++, For[j = 1, j < 5, j++,
       Eins = 
        ReplacePart[
         Eins, {i, j} -> 
          Normal[Normal[
            Series[Ricci[[i, j]] - 1/2 met[[i, j]] RScal, {n,0,1}]]]]]];