The `integral2` function of the **pracma** package is a possibility:

    library(pracma)
    
    f <- function(x ,y) x + y
    
    integral2(f, xmin = 0, xmax = 3, ymin = function(x) x, ymax = function(x) x+2)

___

If you are not allowed to use a package, you can nest the `integrate` function:

    inner <- function(x) integrate(function(y) x+y, lower = x, upper = x + 2)$value
    integrate(Vectorize(inner), lower = 0, upper = 3)

Both methods give `24`, an approximate value of the double integral.

___

Finally, you can use the **SimplicialCubature** package after noticing that the region of integration can be split into two triangles (= simplicies). Moreover, the integrand is a polynomial function, and then **SimplicialCubature** offers the possibility to get the *exact* value of the integral.

    library(SimplicialCubature)
    
    S1 <- cbind(c(0,0), c(3,3), c(0,2)) # first triangle
    S2 <- cbind(c(0,2), c(3,3), c(3,5)) # second triangle
    S <- array(c(S1, S2), dim = c(2, 3, 2))
    
    P <- definePoly(coef = c(1,1), k = cbind(c(1,0), c(0,1)))
    printPoly(P) # x + y
    
    integrateSimplexPolynomial(P, S)
    # 24