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)