x <- seq(0,1,leng=250)
alpha <- 13
beta <- 3

plot(x,dbeta(x,alpha,beta), type="l")
abline(v=qbeta(c(.025,.975),alpha,beta),col="red")

temp <- uniroot( function(x) dbeta(x,alpha,beta) - dbeta( qbeta( .95 +
pbeta(x,alpha,beta), alpha,beta), alpha,beta),
c(0,qbeta(.05,alpha,beta) ))

abline(v=c( temp$root, qbeta(.95+pbeta(temp$root, alpha,beta),
alpha,beta)), col="blue")
abline(h=dbeta(temp$root, alpha,beta), col="green")
cat("Area 1 =",pbeta(temp$root,alpha,beta),"Area 2 =",1-pbeta(qbeta(.95+pbeta(temp$root, alpha,beta),alpha,beta),alpha,beta),"\n")