Solving triple integral with variable bounds based on outer integral

Question

The current code I have is below. The equation I'm trying to solve is below as well. It can be solved in wolfram alpha but I'm unable to get R to perform it.

library(pracma)
fun <- function(x,y) {exp(-(x+y))}
fun2 <- function(z) {integral2(fun, 0, z, 0, 2-z)$Q}
integral(fun2, 0, 2)

Copy this into wolfram alpha to see the pretty text:

integral_0^(1) integral_0^z integral_0^(2-z) e^(-(x+y)) dx dy dz

Answer 1

You were almost there: integral needs a vectorized function.

library(pracma)
fun <- function(x,y) {exp(-(x+y))}
fun2 <- function(z) {integral2(fun, 0, 2-z, 0, z)$Q}
Fun2 <- Vectorize(fun2)
integral(Fun2, 0, 1)
# [1] 0.2706706

Of course, Wolfram Alpha returns a symbolic solution and will thus be more exact. Still, R's result is quite good.