#First, install the package "cubature" with install.packages("cubature").
#You only need to do it once.

library(cubature) # load the package "cubature"
f <- function(x) { 6 * x[1] * (x[2] > x[1])} # This is the density, "x" is a 2D-vector

E.XY <- adaptIntegrate(function(x) x[1]*x[2]*f(x), lowerLimit = c(0, 0), 
                       upperLimit = c(1, 1))
E.X <- adaptIntegrate(function(x) x[1]*f(x), lowerLimit = c(0, 0), 
                      upperLimit = c(1, 1))
E.Y <- adaptIntegrate(function(x) x[2]*f(x), lowerLimit = c(0, 0), 
                      upperLimit = c(1, 1))

cov.XY <- E.XY$integral - E.X$integral * E.Y$integral #E.X, E.Y, E.XY have several fields like 
                                                      #E.X$integral E.X$error and so on. 
                                                      #We need only the integral field.
cov.XY