[R] Improving result from integrate

[. .]

Hi all,

I am utilizing integrate() to get the area below a function that shape
like an exponential, getting very low for higher values of the
horizontal coordinate y.

In the code below, I present the tricky way I can improve the result
obtained. I compare this result with the one I get from Mathematica
and also with the exact solution I have for this particular case.

/*----------
func <- function(y, a, rate){
  x <- function(n,rate) {
    rate*exp(-n*rate)
  }
  boi <- function(y,n,a){
                      w <- y*log(a*n)-lfactorial(y)-a*n
                      exp(w)
                    }
  f <- function(n){
    boi(y,n,a)*x(n,rate)
  }
  r <- 0
                    r1 <- 1
                    x1 <- 0
                    dx <- 20
                    while(r1 > 10e-1000){
                      r1 <- integrate(f,x1,x1+dx)$value
                      r <- r + r1
                      x1 <- x1 + dx
                    }
                    r + integrate(f,x1,Inf)$valu
}
func(200,0.1,0.1)
----------*/

Altought I get better results, the value of dx must be carefully
selected. So I ask, there is another method that can give me better
results?