library(hdrm)

out <- Ex3.1()   # Takes a while, reduce N for quicker results
Fig3.5(out)

Fig3.6(parlist=list(mfrow=c(1,1)))