How can I create a multivariate random normal distribution from separate normal distributions?

Question

In R I can sample from a multivariate normal distribution as follows:

n         <- 100 # sample size
mu        <- c(1, 5) # mean for each distribution           
Va        <- c(0.2, 1) # variation for each distribution        
CORa      <- 0 # correlation between the two                     
Cova      <- CORa*(sqrt(Va[1])*Va[2])  # covariance between the two
VCVmatrix <- matrix(c(Va[1], Cova, Cova, Va[2]), 2, 2) # variance-covariance matrix
MVNData   <- mvrnorm(n, mu, VCVmatrix) # function for multivariate randome normal sampling     

But I want to be able to replicate this without using the mvrnorm function (I've to create these data in a different language that doesn't have this function).

Someone at stackoverflow suggested this but it's not working out exactly right:

rho <- 0
dim1 <- rnorm(100, 1, 0.2)
dim2 <- rho * dim1 + sqrt(1 - rho ^ 2) * rnorm(100, 5, 1)

If I compare the histograms of the distributions from the two methods, dim1 seems to have a much smaller variation than MVNData[,1].

hist(MVNData[,1], col='blue')
hist(dim1, col='red',add=T)

hist(MVNData[,2], col='blue')
hist(dim2, col='red',add=T)

Hope you can help.