Please try:
# bivariate normal distribution plot
library(MASS)
Sigma <- matrix(c(1, 0.5, 0.5, 1), 2)
bivnorm <- mvrnorm(n=1000, rep(0, 2), Sigma)
bivnorm.kde <- kde2d(bivnorm[,1], bivnorm[,2], n=50)
persp(bivnorm.kde, phi=45, theta=40, shade=0.1, main="bivariate normal distribution plot")
image(bivnorm.kde, main="contour plot for bivariate normal distribution")
contour(bivn.kde, add=TRUE)
# end

