\begin{Schunk}
\begin{Sinput}
 x1 <- x2 <- seq(-10, 10, length = 51) 
 dens <- matrix(dmvnorm(expand.grid(x1, x2), sigma = rbind(c(3, 2), c(2, 3))), ncol = length(x1)) 
 s3d <- scatterplot3d(x1, x2, seq(min(dens), max(dens), length = length(x1)), type = "n", grid = FALSE, angle = 70, zlab = expression(f(x[1], x[2])), xlab = expression(x[1]), ylab = expression(x[2]), main = "Bivariate normal distribution")
 for(i in length(x1):1){
     s3d$points3d(rep(x1[i], length(x2)), x2, dens[i,], type = "l")
  }
 for(i in length(x2):1) {
    s3d$points3d(x1, rep(x2[i], length(x1)), dens[,i], type = "l")
  }
\end{Sinput}
\end{Schunk}
