a <- 1 f1 <- function(t) 6*a*exp(-3*a*t) - 6*a*exp(-6*a*t) f2 <- function(t) 12*a*exp(-3*a*t) - 8*a*exp(-4*a*t) - 10*a*exp(-5*a*t) + 6*a*exp(-6*a*t) f3 <- function(t) 24*a*exp(-3*a*t) - 48*a*exp(-4*a*t) + 30*a*exp(-5*a*t) - 6*a*exp(-6*a*t) t <- (0:2000)/1000 plot(t, sapply(t, f1), type="l", xlab="t", ylab="dens(T=t)") text(0.6, 1.5, labels="dens(T1=t)") points(t, sapply(t, f2), type="l") text(0.8, 1.1, labels="dens(T2=t)") points(t, sapply(t, f3), type="l") text(1.4, 0.4, labels="dens(T3=t)")