Gompertz-Makeham Model Life Expectancy using Numerical Integration and Analytical Solution of Castellares et al. (2020) I am trying to calculate life expectancies for the Gompertz-Makeham model, but can't replicate the results of the paper which gives the formulas. I use R to calculate the formula using numerical integration. I provide the formulas and my code below and would be happy if someone could check if I am missing something.

Castellares et al. (2020) provide an analytical solution for the life expectancy of the Gompertz-Makeham model:
$$ \mathcal{M}_{GM} =  \frac{1}{b} \exp\left( \frac{a e^{bx}}{b} \right ) \left( \frac{a e^{bx}}{b} \right)^{c/b} \Gamma\left( - \frac{c}{b}, \frac{a e^{bx}}{b} \right) \qquad \text{with } $$
$$\Gamma(\cdot, \cdot) = \int_x^\infty t^{u-1} e^{-t} dt, \; x>0, \quad u \in \mathbb{R}$$
$\Gamma(\cdot, \cdot)$ is called the upper incomplete gamma distribution. Alternatively, the same result can be achieved via numerical integration of:
$$ \int_0^\infty \exp\left[-ct -   \frac{a e^{bx}}{b} (e^{xb} - 1) \right]
$$
They estimate the parameters $a, b$ and $c$ and denote the age by $x$:

With my implementation of these formulas in R i can't replicate these findings. The analytical solution gives:
# Analytical Solution for age 30
a = 0.00014
b = 0.11521
c = 0.00033
t = 30


gamma_u = -(c/b)              # First parameter for the upper incomplete gamma distribution
gamma_x = (a * exp(b*t)) / b  # Second parameter for the upper incomplete gamma distribution

# Function to integrate over to get the upper incomplete gamma distribution
intfun = function(t, u){
  t^(u-1) * exp(-t)
}

# Result for the upper incomplete gamma distribution
gammainc = integrate(f = intfun,
                     lower = gamma_x,
                     upper = Inf,
                     u = gamma_u)$value

# Analytic formula:
(1/b) * exp(gamma_x) * (gamma_x^(c/b)) * gammainc

# Here the result of the computation:
[1] 24.4011


and the numerical integration gives the same:
## Numerical integration
intfun2 = function(t, a, b, c, x){

  exp( -c*t - (a* exp(b*x))/(b) * (exp(b*t) - 1) )

}

integrate(intfun2,
          lower = 0,
          upper = Inf,
          a = 0.00014, b = 0.11521, c = 0.00033,
          x = 30)

# Here the result of the computation:
24.4011 with absolute error < 5.7e-06


 A: After some research I came upon a GitHub repository of one of the authors in which he used Bell Regression to model mortality. Albeit the parameters $a, b$ and $c$ are derived by different methods as in Castellares et al. (2020), the script also contains code that implements life expectancies using the above formula.
# Makeham
# Gamma-Imcompleta superior
Gamma_Inc_sup = function(s, z){
  integrate(function(t){
    t^(s-1)*exp(-t)
  }, z, Inf)$value
}

ex_mak = function(t, a, b, c){
  (exp(a*exp(b*t)/b)/b)*((a*exp(b*t)/b)^(c/b))*Gamma_Inc_sup(-c/b,a*exp(b*t)/b)
}

This code is equivalent to my code in the above question.
What I found is that the parameter estimation uses mortality data for the ages 30-60, which are the ages in which it is log-linear. Consequently, the code calculates the life expectancy at age  30 using something like ex_mak(0, a, b, c), i.e. t = 0corresponds to age 30.
Using then the values 0, 30 and 60 for the ages 30, 60 and 90 and the parameters as in my above question I find the results:
ex_mak(0, 0.00014, 0.1152, 0.00033)
[1] 52.84945

ex_mak(30, 0.00014, 0.1152, 0.00033)
[1] 24.40485

ex_mak(60, 0.00014, 0.1152, 0.00033)
[1] 4.506404

These are still slightly different to the values in the paper:

but I strongly suppose can attribute these to rounding errors.
