How to obtain the probability distribution of a sum of dependent discrete random variables more efficiently

Question

I hope you are well. I was wondering if you could help me with the question provided in the attached link, please. Below the link I attach an R-code that solves the problem recursively for particular values of the parameters of the distributions involved. However, I realized that this method is inefficient. Thanks a lot for your help.

How to obtain the probability distribution of a sum of dependent discrete random variables more efficiently

library(boot)     # The library boot is necessary to use the command inv.logit.

TMax <- 500       # In this R-code, I am using TMax instead of using T.
M <- 2000
beta0 <- 1
beta1 <- 0.5 
Prob_S <- function(k, r){        # In this R-code, I am using r instead of using t.
    if(r == 1){
        Aux <- dbinom(x = k, size = M, prob = inv.logit(beta0))
        }
    if(r %in% 2:TMax){
        Aux <- 0
        for(u in 0:k){
            Aux <- Aux + dbinom(x = k - u, size = M - u, 
                prob = inv.logit(beta0 + beta1 * u)) * Prob_S(u, r - 1)
            }
        }
    Aux
    }

m <- 300
P <- Prob_S(k = m, r = TMax)    # Computing P takes a loooong time.   :(