diff --git a/R/ByChunks.R b/R/ByChunks.R
index 37a554c1e44d8fda4bc60a6a2a0732c54e4009b3..25876082efc960f20f1b08b453d5b8e47a434119 100644
--- a/R/ByChunks.R
+++ b/R/ByChunks.R
@@ -735,6 +735,9 @@ ByChunks <- function(step_fun, cube_headers, ..., chunks = 'auto',
       }
       #TODO: Find a better way to assign the names of data. When multiple steps for Compute is available, this way may fail.
       names(data) <- names(cube_headers)
+      # Load libraries
+      invisible(lapply(attr(step_fun, 'UseLibraries'), library, character.only = TRUE))
+
       t_begin_compute <- Sys.time()
       result <- multiApply::Apply(data, 
                                   target_dims = attr(step_fun, 'TargetDims'), 
diff --git a/tests/testthat/test-Compute-timedim.R b/tests/testthat/test-Compute-timedim.R
index 80d96ff0bf2ade58ae1b30c9f2cadd1c6e6e8e36..6b8138963f35099da3a455a4d08793bd30917827 100644
--- a/tests/testthat/test-Compute-timedim.R
+++ b/tests/testthat/test-Compute-timedim.R
@@ -16,21 +16,37 @@ suppressWarnings(
                 retrieve = FALSE)
 )
   
+#------------(1) specify package name in front of function
   fun_spring <- function(x) {
-      y <- s2dv::Season(x,  time_dim = 'time', monini = 1, moninf = 3, monsup = 5)
-      return(y)
+    y <- s2dv::Season(x,  time_dim = 'time', monini = 1, moninf = 3, monsup = 5)
+    dummy <- pryr::object_size(y)
+    return(y)
   }
-  
   step1 <- Step(fun = fun_spring,
                target_dims = c('var', 'time'),
                output_dims = c('var', 'time'))
-  
   wf1 <- AddStep(data, step1)
 
+#------------(2) use use_libraries 
+  fun_spring2 <- function(x) {
+    y <- Season(x,  time_dim = 'time', monini = 1, moninf = 3, monsup = 5)
+    dummy <- object_size(y)
+    return(y)
+  }
+  step2 <- Step(fun = fun_spring2,
+               target_dims = c('var', 'time'),
+               output_dims = c('var', 'time'),
+               use_libraries = c('s2dv', 'pryr'))
+  wf2 <- AddStep(data, step2)
+
+#-------------------------------------
 suppressWarnings(
-  res1 <- Compute(wf1,
-              chunks = list(ensemble = 2))
+  res1 <- Compute(wf1, chunks = list(ensemble = 2))
 )
+suppressWarnings(
+  res2 <- Compute(wf2, chunks = list(ensemble = 2))
+)
+#=======================================
 
 expect_equal(
 attr(data, 'Dimensions'),
@@ -51,7 +67,10 @@ c(256.4469, 260.3636),
 tolerance = 0.0001
 )
 
-
+expect_equal(
+c(res1$output),
+c(res2$output)
+)
 
 
 })