I create a function to output the mean of samples

mean.dr <- function (n, m){
  mtx <- matrix (nrow = n, ncol = m)
  
  mtx[1:n,] <- replicate (m, rchisq (n = n, df = 3))
  
  M <- colMeans (mtx)
  
  meanM <- mean(M)
  sdM <- sd(M)
  
  
  # save them as "result"
  result_toplot <- as.numeric (M)
  
  # plot these means as a histogram
  hist (result_toplot, probability = TRUE,
        col = "cornflowerblue", border = "darkgreen",
        main = paste("histogram of means ( n =", n,   ", df = "   ,3,   ")"))
  
  x <- seq (min(result_toplot), max(result_toplot), 
            length.out = 1000)
  mean_result <- meanM
  sd_result <- sqrt (6 / n)
  
  lines (x, dnorm (x, mean = mean_result, sd = sd_result),
         col = "red", lwd = 3)
}

Try it:

mean.dr (n = 500, m = 1000)

Use the function four times

par (mfrow = c (2, 2))

mean.dr (n = 100, m = 1000)
mean.dr (n = 200, m = 1000)
mean.dr (n = 100, m = 2000)
mean.dr (n = 200, m = 2000)

par (mfrow = c (1, 1))

You can also use Loop function:

par (mfrow = c (2, 2))

ns <- c (100, 200)
ps <- c (1000, 2000)

for (n in ns) {
  for (p in ps) {
    mean.dr(n, p)
  }
}

par (mfrow = c (1, 1))