Does this code demonstrate the central limit theorem? Does this code demonstrate the central limit theorem? This is not a homework assignment! Au contraire, I'm a faculty teaching some methods to non-stats students. 
library(tidyverse)
#Make fake data
population<-rnorm(1000000, mean=100, sd=10)

#Draw 100 samples of size 5
map(1:100, ~sample(population, size=5)) %>% 
  #calculate their mean
  map(., mean) %>% 
  #unlist 
  unlist() %>% 
  #draw histogram of sample means
  hist(, xlim=c(80,120))

#Repeat but with sample size 500
map(1:100, ~sample(population, size=500)) %>% 
  map(., mean) %>% 
  unlist() %>% 
  hist(., xlim=c(80,120))

#Repeat but with sample size 1000
map(1:100, ~sample(population, size=1000)) %>% 
  map(., mean) %>% 
  unlist() %>% 
  hist(., xlim=c(80,120))




 A: Here's an example of one of my suggestions from comments. Means of samples of size n=100000 (takes about 20 seconds or so, be patient):
  ln.mean = replicate(1000,mean(rlnorm(100000,0,4)))
  hist(ln.mean,n=100)


Even at this huge sample size, the distribution of sample means is still really skew -- but the central limit theorem nevertheless applies here - even the "classic" CLT.
A: Maybe use something like the following (simpler, more direct) R code to show that
averages of a dozen standard uniform random variables are
difficult to distinguish from normal.
set.seed(1126)
a = replicate(5000, mean(runif(12))
shapiro.test(a)

        Shapiro-Wilk normality test

data:  a
W = 0.99965, p-value = 0.565

plot(qqnorm(a))


Then use R code to show that averages of 50, or even 100, standard exponential
random variables are easy to distinguish from normal. What is the distribution
of $A = \bar X_{100}?$
set.seed(1127)
a = replicate(5000, mean(rexp(100)))
shapiro.test(a)$p.val
 [1] 1.675877e-06

However, averages of 1000 standard exponentials are more difficult to distinguish from normal.
set.seed(1127)
a = replicate(5000, mean(rexp(1000)))
shapiro.test(a)$p.val
[1] 0.2413559

A: Here's a complete study in a few lines. 
For a given set of sample sizes n and underlying distribution r, it generates n.sim independent samples of each size from that distribution, standardizes the empirical distribution of their means, plots the histogram, and overplots the standard Normal density in red.  The CLT says that when the underlying distribution has finite variance, the red curve more and more closely approximates the histogram.

The first three rows illustrate the process for sample sizes of $10,20,100,500$ and underlying Normal, Gamma, and Bernoulli distributions.  As sample size increases the approximation grows noticeably better.  The bottom row uses a Cauchy distribution.  Because a key assumption of the CLT (finite variance) does not hold in this case, its conclusion doesn't hold, which is pretty clear.
Execution time is about one second.
f <- function(n, r=rnorm,  n.sim=1e3, name="Normal", ...) {
  sapply(n, function(n) {
    x <- scale(colMeans(matrix(r(n*n.sim, ...), n))) # Sample, take mean, standardize
    hist(x, sub=name, main=n, freq=FALSE, breaks=30) # Plot distribution
    curve(dnorm(x), col="Red", lwd=2, add=TRUE)      # Compare to standard Normal
  })
}
n <- c(5,20,100,500)
mfrow.old <- par(mfrow=c(4,length(n)))
f(n)
f(n, rgamma, shape=1/2, name="Gamma(1/2)")
f(n, function(n) runif(n) < 0.9, name="Bernoulli(9/10)")
f(n, rt, df=1, name="Cauchy")
par(mfrow=mfrow.old)

