MCMC bivariate Normal

"MCMC" <-
function(my,sy,mx,sx,r,n,x0=mx,pts=TRUE)
{
# MCMC generation of bivariate normals
X <- rep(0,n+1)
Y <- rep(0,n+1)
X[1] <- x0
sr <- sqrt(1-r^2)
if(pts) plot(X,Y,xlim=c(mx-4*sx,mx+4*sx),ylim=c(my-4*sy,my+4*sy),type="n")
for(i in 1:n)
{
      Y[i] <- rnorm(1,my + (r*sy/sx)*(X[i]-mx),sy*sr)
      X[i+1] <- rnorm(1,mx + (r*sx/sy)*(Y[i]-my),sx*sr)
      if(pts) points(X[i],Y[i],pch=20,col=rgb(i/n,0,0))
      }
cbind(X,Y)
}