# you can go to https://cran.r-project.org/
# and use the search facility to look for
# R packages other people have already written

# e.g., searching on 'multivariate normal'
# gets you the package mvtnorm

# first you need to install a package;
# then you can load it with a command of the form

library( mvtnorm )

n.grid <- 100

x.1.grid <- seq( -3, 3, length = n.grid )

x.2.grid <- seq( -3, 3, length = n.grid )

mu <- c( 0, 0 )

print( Sigma <- matrix( c( 1, 0.9, 0.9, 1 ), 2, 2 ) )

#      [,1] [,2]
# [1,]  1.0  0.9
# [2,]  0.9  1.0

bivariate.normal.density <- matrix( NA, n.grid, n.grid )

for ( i in 1:n.grid ) {

  for ( j in 1:n.grid ) {

    bivariate.normal.density[ i, j ] <- dmvnorm( c( x.1.grid[ i ], 
      x.2.grid[ j ] ), mu, Sigma )

  }

}

persp( x.1.grid, x.2.grid, bivariate.normal.density,
  xlab = 'x1', ylab = 'x2', zlab = 'Density' )

contour( x.1.grid, x.2.grid, bivariate.normal.density,
  xlab = 'x1', ylab = 'x2' )

par( mfrow = c( 2, 2 ) )

persp( x.1.grid, x.2.grid, bivariate.normal.density,
  xlab = 'x1', ylab = 'x2', zlab = 'Density' )

persp( x.1.grid, x.2.grid, bivariate.normal.density,
  xlab = 'x1', ylab = 'x2', zlab = 'Density',
  theta = 0, phi = 60 )

persp( x.1.grid, x.2.grid, bivariate.normal.density,
  xlab = 'x1', ylab = 'x2', zlab = 'Density',
  theta = 30, phi = 15 )

persp( x.1.grid, x.2.grid, bivariate.normal.density,
  xlab = 'x1', ylab = 'x2', zlab = 'Density',
  theta = 30, phi = 60 )

par( mfrow = c( 1, 1 ) )

library( plot3D )

persp3D( x.1.grid, x.2.grid, bivariate.normal.density,
  xlab = 'x1', ylab = 'x2', zlab = 'Density' )