R code used to produce the graph shown below.

library(gplots) 
data(airquality) 
# Compute means and standard deviations for each month 
mean.values <- by(airquality$Temp, airquality$Month, mean) 
sd.values <- by(airquality$Temp, airquality$Month, sd) 
barplot2(mean.values, ylim=c(0,99), 
         xlab="Month", ylab="Temperature (F)", 
         plot.ci=TRUE, 
         ci.u = mean.values+1.96*sd.values, 
         ci.l = mean.values-1.96*sd.values, ci.lwd=3)