calculate the gradient for the multi-variate normal distribution

multinorm <- function(x, y, mean = c(0, 0), cov = matrix(c(1, 0, 0, 1), nrow = 2)) 
{
  x <- as.vector(x)
  y <- as.vector(y)
  mean <- as.vector(mean)
  cov <- as.matrix(cov)
  
  n <- length(x)
  if (n != length(y)) stop("data vectors must have the same length")
  if (n != length(mean)) stop("mean vector must have the same length as data vectors")
  if (nrow(cov) != n) stop("covariance matrix must have the same number of rows as the length of mean and data vectors")
  if (nrow(cov) != ncol(cov)) stop("covariance matrix must be square")
  
  exx <- cov[1, 1]
  eyy <- cov[2, 2]
  exy <- cov[1, 2]
  e <- 2 * exy
  d <- exx