diff --git a/R/expected_event.R b/R/expected_event.R
index c105678f..efe68079 100644
--- a/R/expected_event.R
+++ b/R/expected_event.R
@@ -165,22 +165,13 @@ expected_event <- function(
 
   # Create 3 step functions (sf) ----
   # Step function to define enrollment rates over time
-  sf_enroll_rate <- stats::stepfun(c(0, cumsum(enroll_rate$duration)),
-    c(0, enroll_rate$rate, 0),
-    right = FALSE
-  )
+  sf_enroll_rate <- stepfun2(c(0, cumsum(enroll_rate$duration)), c(0, enroll_rate$rate, 0))
   # step function to define failure rates over time
   start_fail <- c(0, cumsum(fail_rate$duration))
   fail_rate_last <- nrow(fail_rate)
-  sf_fail_rate <- stats::stepfun(start_fail,
-    c(0, fail_rate$fail_rate, fail_rate$fail_rate[fail_rate_last]),
-    right = FALSE
-  )
+  sf_fail_rate <- stepfun2(start_fail, c(0, fail_rate$fail_rate, fail_rate$fail_rate[fail_rate_last]))
   # step function to define dropout rates over time
-  sf_dropout_rate <- stats::stepfun(start_fail,
-    c(0, fail_rate$dropout_rate, fail_rate$dropout_rate[fail_rate_last]),
-    right = FALSE
-  )
+  sf_dropout_rate <- stepfun2(start_fail, c(0, fail_rate$dropout_rate, fail_rate$dropout_rate[fail_rate_last]))
 
   # combine sub-intervals from enroll + failure + dropout #
   # impute the NA by step functions
@@ -223,7 +214,7 @@ expected_event <- function(
   if (simple) {
     ans <- sum(df$nbar)
   } else {
-    sf_start_fail <- stats::stepfun(start_fail, c(0, start_fail), right = FALSE)
+    sf_start_fail <- stepfun2(start_fail, c(0, start_fail))
     ans <- data.frame(
       fail_rate = df$fail_rate_var,
       event = df$nbar,
diff --git a/R/utils.R b/R/utils.R
index fad1d1c5..40d0f4bb 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -48,12 +48,21 @@ is_wholenumber <- function (x, tol = .Machine$double.eps^0.5)  {
   abs(x - round(x)) < tol
 }
 
-# a faster version of stats::stepfun() since we don't need to consider interpolation
+# a faster version of stats::stepfun() since we don't need to consider interpolation;
+# right = FALSE means left-closed intervals [x0[i], x0[i+1]), i.e., the function jumps
+# at x0[i] (right-continuous); right = TRUE means right-closed (x0[i], x0[i+1]]
 stepfun2 <- function(x0, y, right = FALSE) {
   x0; y  # avoid lazy evaluation: evaluate right now
-  function(x) {
-    i <- findInterval(x, x0, left.open = right)
-    y[i + 1]
+  # for small breakpoint vectors, a for-loop of vectorized comparisons (x >= b)
+  # is faster than findInterval(), which has dispatch overhead and uses binary
+  # search (overkill for 3-5 breakpoints)
+  if (length(x0) <= 10L) {
+    if (right)
+      function(x) { i <- 1L; for (b in x0) i <- i + (x > b); y[i] }
+    else
+      function(x) { i <- 1L; for (b in x0) i <- i + (x >= b); y[i] }
+  } else {
+    function(x) y[findInterval(x, x0, left.open = right) + 1L]
   }
 }