R-coding HPD for the beta-binomial model

r1 <- qbeta(0.025,3,9)
r2 <- qbeta(0.975,3,9)
r2-r1

f <- function(x,k) dbeta(x,3,9)-k
ktest <- seq(0.4,0.8,0.01)
res <- rep(0,length(ktest))
l <- 0
for (i in ktest)
{
  l <- l+1
  r1 <- uniroot(f,c(0,0.2),k=i)$root
  r2 <- uniroot(f,c(0.2,0.8),k=i)$root
  res[l] <- pbeta(r2,3,9)-pbeta(r1,3,9)
}

res[order(abs(res-0.95))[1]]
kalpha <- ktest[order(abs(res-0.95))[1]]

r1 <- uniroot(f,c(0,0.2),k=kalpha)$root
r2 <- uniroot(f,c(0.2,0.8),k=kalpha)$root
r2-r1