set.seed(12345)

n<-10000;
z<-rnorm(n,0,1)
par(mfrow=c(1,2))

hist(z,  freq=F, col="grey", xlab="", xlim=c(-4, 4), breaks=100)
par(new=T)
plot(density(z), axes=F, main="", xlim=c(-4, 4), lwd=2, col="blue")
curve(dnorm(x,0,1),xlim=c(-4, 4), lwd=2)