m <- matrix(1:15, 3, 5)
cn <- combn(1:ncol(m), 2)
res <- apply(cn, 2, function(x) {
  m[, x[1]] * m[, x[2]]
})
colnames(res) <- paste(cn[1, ], cn[2, ], sep = "_")
print(m)
print(res)