myfun5 <-  function(x, theta, epi) {
if(epi <= 0) stop("epi must > 0")
u <- (x-theta)*sqrt(5)/epi
INDE <- (u < sqrt(5)) & (u > -sqrt(5))
u[u >= sqrt(5)] <- 0
u[u <= -sqrt(5)] <- 1
y <- 0.5 - (u - (u)^3/15)*3/(4*sqrt(5))
u[ INDE ] <- y[ INDE ]
return(u)
}