Here are 4 quick ways to plot 2D distributions.  See here for more color options.

library(MASS)   #kde2d() and mvrnorm()
VV=matrix(c(10,8,8,10),2,2)
coef=mvrnorm(10000,c(2,10),VV)
par(mfrow=c(2,2), mar=c(4,4,1,1))
### 1
plot(coef[,1],coef[,2])
### 2
contour(kde2d(coef[,1],coef[,2])) 
### 3
persp(kde2d(coef[,1],coef[,2])) 
### 4 
image(kde2d(coef[,1],coef[,2], n = 100))
box()

This code adds estimated bands to 2D plots.  This can be used if you have posterior draws for two parameters (ie. from MCMC).  Estimated confidence bands of any level can be specified; the default is 68% and 95%.  This function takes some time; if you have a lot of points you may want to thin them before plotting (ie. use every 25th).

This is a plot I made and used quite a bit.  The input is a matrix of residuals (subject x number of MCMC iterations) coming from MCMC.  Here is how I made this plot.