a <- 1
F1 <- function(t) 1 - 2*exp(-3*a*t) + exp(-6*a*t)
F2 <- function(t) 1 - 8*exp(-3*a*t) + 12*exp(-4*a*t) - 6*exp(-5*a*t) + exp(-6*a*t)
f1 <- function(t) 6*a*exp(-3*a*t) - 6*a*exp(-6*a*t)
f2 <- 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, f1(t), type="l", xlab="t", ylab="dens(T=t)")
points(t, f2(t), type="l")