Random effects:
Groups   Name Variance  Std.Dev. Corr
subject  X21  8.558e+00 2.925380
         X22  2.117e-03 0.046011 -1.00
         X23  2.532e-05 0.005032  1.00 -1.00
Residual      1.453e+00 1.205402
Number of obs: 100, groups:  subject, 20

> Omega
          [,1] [,2] [,3]
[1,] 0.6181442    0    0
[2,] 0.0000000    0    0
[3,] 0.0000000    0    0

M1 <- lmer(ym ~ 0+XC1 + (0+X2 | subject))
rr <- ranef(M1,condVar=TRUE)

pv <- attr(rr[[1]],"postVar")
str(pv)

> pv[,,1]
              [,1]          [,2]          [,3]
[1,]  0.2913922395 -4.588735e-03  5.017337e-04
[2,] -0.0045887347  7.523883e-05 -8.101042e-06
[3,]  0.0005017337 -8.101042e-06  8.773362e-07

N=20          #number of subject
n=5           #number of observations per subject

sigma2 =1

#equally space point not include the start point
time <- function(from, to, length.out) {
    length.out <- length.out + 1
    result <- seq(from, to, length.out = length.out)
    result <- result[-1]
    return(result)
}

subject = matrix(0,nrow=N*n,ncol=1)
for(i in 1:N){
    for(j in (n*(i-1)+1):(n*i)){
        subject[j]=i
    }
}


X = array(0, dim = c(N, n, 3))#each X[i,,] is a nx3 matrix
for (i in 1:N){
    for (j in 1:n){
    X[i,j,] <-c(1,time(0,10,n)[j],(time(0,10,n)[j])^2)
    }
}


y = array(0, dim = c(N, n, 1))


Omega <- matrix(0,nrow=3,ncol=3)
Omega[1,1] = runif(1,0.01,1.01)#only omega1^2 is not equal to 0

beta <-rep(0,5)
beta[1]= rnorm(1,mean=0.01,sd=1) #mu0
beta[2]= rnorm(1,mean=0.005,sd=1) #mu1
beta[3]= rnorm(1,mean=0.0025,sd=1) #mu2


C1 = array(0, dim = c(N, 3, 5))
for(i in 1:N){
    C1[i,1,1]=C1[i,2,2]=C1[i,3,3]=1
}


muy = array(0, dim = c(N, n, 1))     #store the expextation of each yi
Cov = array(0, dim = c(N, n, n))     #store the covariance matrix of y

for (i in 1:N){
    muy[i,,] <- X[i,,]%*%C1[i,,]%*%beta
    Cov[i,,] <- X[i,,]%*%Omega%*%t(X[i,,])+ sigma2*diag(n)
    y[i,,] <- mvrnorm(n = 1, muy[i,,], Cov[i,,])
}

ym <- as.vector(y[,,1])


#change X into X2, which is in a matrix format, easy for compitation later
X2 <- rbind(X[1,,],X[2,,])
for(i in 2:(N-1)){
    X2 = rbind(X2,X[i+1,,])
}


XC1=matrix(0,nrow=N*n,ncol=5)
for(i in 1:N){
XC1[(n*(i-1)+1):(i*n),]=X[i,,]%*%C1[i,,]
}



M1<-lmer(ym ~ 0+XC1+(0+X2|subject))