$$ \begin{align} \alpha_m^l &= \sum_{m=0}^{l-1}{\frac{\Gamma(cm+1)\Gamma(cl-cm+1)}{\Gamma(m+1)\Gamma(l-m+1)}} \\ \alpha_j^0 &= \frac{\Gamma(cj+1)}{\Gamma(j+1)}, \quad j=0,1,2,\ldots \\ \alpha_j^{n+1} &= \sum_{m=n}^{j-1}{\alpha_m^n \frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}} \end{align} $$

You can sensibly sum $n$ only up to whatever is an extremely unlikely number of failures within time $T_m$. Choosing increasing values of $j$ will result in progressively more accurate estimation, so this should be summed up as appropriate for your desired accuracy.

This can all be implemented in Excel (use =EXP(GAMMALN(X)) for $\Gamma(X)$). I would suggest creating $J \times N$ grids (where $J$ and $N$ are your choice for how high to sum over $j$ and $n$) for $\Pr(N(T_m)=n)$ and calculation of $E(N(T_m))$, drawing from vectors created for $\alpha_m^n$ and $\alpha_j^n$.

I'm afraid all the legwork is left for you to do, but hopefully this is a sufficiently comprehensive answer to get you on your way!

$$ \begin{align} \alpha_j^0 &= \frac{\Gamma(cj+1)}{\Gamma(j+1)}, \quad j=0,1,2,\ldots \\ \alpha_j^{n+1} &= \sum_{m=n}^{j-1}{\alpha_m^n \frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}} \end{align} $$

You can sensibly sum $n$ only up to whatever is an extremely unlikely number of failures within time $T_m$. Choosing increasing values of $j$ will result in progressively more accurate estimation, so this should be summed up as appropriate for your desired accuracy.

This can all be implemented in Excel (use =EXP(GAMMALN(X)) for $\Gamma(X)$). I would suggest creating $J \times N$ grids (where $J$ and $N$ are your choice for how high to sum over $j$ and $n$) for $\Pr(N(T_m)=n)$ and calculation of $E(N(T_m))$, drawing from vectors created for $\alpha_m^n$ and $\alpha_j^n$.

I'm afraid all the legwork is left for you to do, but hopefully this is a sufficiently comprehensive answer to get you on your way!

This can all be implemented in Excel (use =EXP(GAMMALN(X)) for $\Gamma(X)$). I would suggest creating $J \times N$ grids (where $J$ and $N$ are your choice for how high to sum over $j$ and $n$) for $\Pr(N(T_m)=n)$ and calculation of $E(N(T_m))$, drawing from vectors created for $\alpha_m^n$ and $\alpha_j^n$.

I'm afraid all the legwork is left for you to do, but hopefully this is a sufficiently comprehensive answer to get you on your way!

EDIT

As stated in a comment below, the notation employed by the authors is a little confusing, here is a clarification with some of my own notation:

1. Define base case $A(0, j) = \frac{\Gamma(cj+1)}{\Gamma(j+1)}$
2. Define recursive case $A(n+1, j) = \sum_{m=n}^{j-1}{A(n, m)\frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}}$
3. Obtain results:
   • $C_n(t) = \Pr(N(t)=n) = \sum_{j=n}^\infty{\frac{(-1)^{j+n}(\lambda t^c)^j A(n, j)}{\Gamma(cj+1)}}$
   • $E[N(t)] = \sum_{n=0}^\infty{n C_n(t)}$

I have not yet coded this as an Excel spreadsheet (though it is doubtless possible), instead here is Python v3 code:

from math import exp
from scipy.special import gamma, gammaln

rate = 1.5
shape = 0.8

Acache = dict()

def A(n, j):
    if not (n, j) in Acache:
        if n == 0:
            Acache[(n, j)] = exp(gammaln(shape*j+1) - gammaln(j+1))
        else:
            Acache[(n, j)] = sum( A(n-1, m) * exp(gammaln(shape*(j-m)+1) - gammaln(j-m+1)) for m in range(n-1, j) )
    return Acache[(n, j)]

def Cn(t, n, J=20):
    return sum( ( (-1)**(j+n) * (rate * t**shape)**j * A(n, j) )/gamma(shape*j+1) for j in range(n, n+J) )

def E(t, J=20, N=20):
    return sum( (n * Cn(t, n, J)) for n in range(0, N) )

C = list(Cn(1, n) for n in range(0, 20))

print('n', 'Cn(1)', sep='\t')
for n in range(0, 20):
    print(n, C[n], sep='\t')
    
print(r'\sum_{n=0}^{19}{C_n(1)} =', sum(C))
    
print('E[N(1)] =', E(1))

Which gives the following output:

n   Cn(1)
0   0.223130160148
1   0.285693698536
2   0.232154412273
3   0.142130392306
4   0.0707688290882
5   0.0299340192311
6   0.0110641801411
7   0.00364513700374
8   0.00108625766843
9   0.000296143808589
10  7.45318952329e-05
11  1.74439294931e-05
12  3.81998309182e-06
13  7.86741813662e-07
14  1.53063575444e-07
15  2.82382329525e-08
16  4.95652671252e-09
17  8.30172205633e-10
18  1.33028122034e-10
19  2.0441696595e-11
\sum_{n=0}^{19}{C_n(1)} = 0.999999999996
E[N(1)] = 1.71339144029