fun1 <- function(x,y) { x^2 +y^2 + 2 } x <- runif(100000,-1000,1000) y <- runif(100000,-1000,1000) result <- fun1(x,y) hist(result, breaks = 100, col=c(blues9))