General ideas are from [here][1] and [here][2].
#R code

1. You first need to simulate a vector of uncorrelated Gaussian random variables, $\bf Z $
        
        NVariables=5
        VariableLen=1000
        Z=matrix(rnorm(NVariables*VariableLen), ncol=NVariables)

2. Create covariance matrix $\Sigma$. Let all variables be correlated with neighbor as 0.5.

        Sigma=matrix(data=0, ncol=NVariables-1,nrow=NVariables-1)
        diag(Sigma)<-0.5
        Sigma=cbind(matrix(data=0,nrow=NVariables-1),Sigma)
        Sigma=rbind(Sigma,matrix(data=0,ncol=NVariables))
        diag(Sigma)<-1
        
3. Find a square root of $\Sigma$. Cholesky decomposition is common choice

        C=chol(Sigma)

4. To obtain random variables with given correlation matrix $\Sigma$ multiply $\bf C$  and $\bf Z$ 

        Y=Z%*%C
5. Use inverse CDF method to obtain any distribution You wish. Here it is lognormal

        Ylog=qlnorm(pnorm(Y))

#Results
[![enter image description here][3]][3]

**Correlation Matrix**


                [,1]        [,2]        [,3]        [,4]        [,5]
    [1,]  1.00000000  0.52817152 -0.01887624 -0.07113405 -0.05551355
    [2,]  0.52817152  1.00000000  0.49392903 -0.03233261 -0.01504632
    [3,] -0.01887624  0.49392903  1.00000000  0.50604908  0.04029076
    [4,] -0.07113405 -0.03233261  0.50604908  1.00000000  0.49000229
    [5,] -0.05551355 -0.01504632  0.04029076  0.49000229  1.00000000


  [1]: https://math.stackexchange.com/questions/446093/generate-correlated-normal-random-variables
  [2]: http://stats.stackexchange.com/questions/184325/how-does-the-inverse-transform-method-work
  [3]: https://i.sstatic.net/OKjTf.png