# A reliability function
a <- 1
F <- function(t) 1 - 3*exp(-a*t) + 3*exp(-2*a*t) - exp(-3*a*t)
f <- function(t) 3*a*exp(-a*t) - 6*a*exp(-2*a*t) + 3*a*exp(-3*a*t)
t <- (0:10000)/1000
plot(t, F(t), type="l", xlab="t", ylab="P(T<=t)")