set.seed(12343)

alpha <- 0.05 
normval <- qnorm(1 - alpha/2) 
numsamp <- 100
numsim <- 10000  
normmat <- matrix(0, nrow = numsim, ncol = 2) 
y <- 1:numsim
ymat <- rbind(y, y) 

for (i in 1:numsim) { 
  samp<-rnorm(numsamp,10,2) 
  sampmean <- mean(samp) 
  sampse <- sqrt(var(samp) / numsamp) 
  normmat[i, ] <- c(sampmean-normval*sampse, sampmean+normval*sampse) 
} 

matplot(t(normmat[1:100,]), ymat[,1:100] , pch = " ", yaxt = "n", ylab = "", xlab="", main="95% Confidennce Interval") # empty plot 
matlines(t(normmat[1:100,]), ymat[,1:100], lty = rep(1, numsim), col = 1) 
abline(v = 10)