windows(5,3)
par(mar=c(3,3,1,1))

#cook up some trajectories like you might get from mcmc
Years = 1980:2000
Nyears = length(Years)
Nvals = 1000
biomat = matrix(rnorm(Nyears*Nvals)+sin(Years/pi)+5, Nvals, Nyears, byrow=T)

CIlevels = seq(.05,.95,.1) #Probability intervals
Nlevels = length(CIlevels)
CIprobs = sort(c(1 - CIlevels/2, CIlevels/2)) #quantiles containing the prob. intervals.

bioquant = apply(biomat,2,quantile, probs=c(.025,.5,.975)) #get median and 95% values
bioquant2 = apply(biomat,2,quantile, probs=CIprobs) #get values for all quantiles
plot(0, type='n', xlim=range(Years), ylim=c(0,1.1*max(bioquant)), xaxs='i', yaxs='i', xlab='', ylab='', main='')
for(ilevel in 1:Nlevels)
{
  polygon(Years[c(1:Nyears,Nyears:1)], c(bioquant2[ilevel,1:Nyears],
    bioquant2[2*Nlevels+1-ilevel,Nyears:1]), border=NA, col=paste('grey',floor(45+50*CIlevels[Nlevels+1-ilevel])))
}
matplot(Years,t(bioquant), type='l', col=1, lwd=c(1,3,1), lty=c(3,1,3), add=T)
box()