par(mfrow = c(2, 2))

norm = density(rnorm(1e6,0,1))
n.minx = min(norm$x)
n.maxx = max(norm$x)
n.miny = min(norm$y)
n.maxy = max(norm$y)

for (i in c(10, 20, 50, 100) ) {
  d = density(rnorm(i,0,1))
  plot ( d, 
         xlim=c( min(n.minx,min(d$x)),max(n.maxx,max(d$x)) ),
         ylim=c( min(n.miny,min(d$y)),max(n.maxy,max(d$y)) )  )
  lines ( norm, col="Red" )
}

Deviations from normal distribution based on number of observations

R