ONES = J(N, 1, 1);
meanvec = (1/N)*t(X)*ONES;
mean_matrix = ONES*t(meanvec);
cov_matrix = (1/n) * t(X- mean_matrix) * (x - mean_matrix);

ONES <- matrix(1, nrow=N, ncol=1)
meanvec <- (1/N) * t(X) %*% ONES
mean_matrix <- ONES %*% t(meanvec)
cov_matrix <- (1/N) * t(X - mean_matrix) %*% (X - mean_matrix)

X
     [,1] [,2] [,3]
[1,]   90   60   90
[2,]   90   90   30
[3,]   60   60   60
[4,]   60   60   90
[5,]   30   30   30

cov_matrix
     [,1] [,2] [,3]
[1,]  504  360  180
[2,]  360  360    0
[3,]  180    0  720

cov(X)
     [,1] [,2] [,3]
[1,]  630  450  225
[2,]  450  450    0
[3,]  225    0  900