ZLHS2016:

# 3-D plots
#
#
mu1<-0 # setting the expected value of x1
mu2<-0 # setting the expected value of x2
s11<-10 # setting the variance of x1
s12<-15 # setting the covariance between x1 and x2
s22<-10 # setting the variance of x2
rho<-0.5 # setting the correlation coefficient between x1 and x2
x1<-seq(-10,10,length=41) # generating the vector series x1
x2<-x1 # copying x1 to x2

f<-function(x1,x2)
{
  term1<-1/(2*pi*sqrt(s11*s22*(1-rho^2)))
  term2<--1/(2*(1-rho^2))
  term3<-(x1-mu1)^2/s11
  term4<-(x2-mu2)^2/s22
  term5<--2*rho*((x1-mu1)*(x2-mu2))/(sqrt(s11)*sqrt(s22))
  term1*exp(term2*(term3+term4-term5))
} # setting up the function of the multivariate normal density
#
z<-outer(x1,x2,f) # calculating the density values
#
persp(x1, x2, z,
      main="Two dimensional Normal Distribution",
#      sub=expression(italic(f)~(bold(x)) ==
#                      frac(1,2~pi~sqrt(sigma[11]~sigma[22]~(1-rho^2))) ~
#                       phantom(0)^bold(.)~exp~bgroup("{",
#                                                     list(-frac(1,2(1-rho^2)),
#                                                      bgroup("[", frac((x[1]~-~mu[1])^2, sigma[11])~-~2~rho~frac(x[1]~-~mu[1],                                                                                                                                                                                 sqrt(sigma[11]))~ frac(x[2]~-~mu[2],sqrt(sigma[22]))~+~
#                                                      frac((x[2]~-~mu[2])^2, sigma[22]),"]")),"}")),
      col="lightgreen",
      theta=30, phi=20,
      r=50,
      d=0.1,
      expand=0.5,
      ltheta=90, lphi=180,
      shade=0.75,
      ticktype="detailed",
      nticks=5) # produces the 3-D plot
#
mtext(expression(list(mu[1]==0,
                      mu[2]==0,
                      sigma[11]==10,
                      sigma[22]==10,
                      sigma[12]==15,
                      rho==0.5)),
      side=3) # adding a text line to the graph