diff --git a/DESCRIPTION b/DESCRIPTION
index f26ac97d..4ff62243 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: PKPDsim
 Type: Package
 Title: Tools for Performing Pharmacokinetic-Pharmacodynamic Simulations
-Version: 1.4.1
-Date: 2025-04-09
+Version: 1.5.0
+Date: 2025-07-08
 Authors@R: c(
         person("Ron", "Keizer", email = "ron@insight-rx.com", role = c("aut", "cre")),
         person("Jasmine", "Hughes", email = "jasmine@insight-rx.com", role = "aut"),
diff --git a/NAMESPACE b/NAMESPACE
index 145be003..92cbb2de 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -15,7 +15,6 @@ export(advan_create_data)
 export(advan_parse_output)
 export(advan_process_infusion_doses)
 export(apply_duration_scale)
-export(apply_lagtime)
 export(available_default_literature_models)
 export(calc_auc_analytic)
 export(calc_ss_analytic)
diff --git a/R/apply_lagtime.R b/R/apply_lagtime.R
deleted file mode 100644
index 85718ccf..00000000
--- a/R/apply_lagtime.R
+++ /dev/null
@@ -1,40 +0,0 @@
-#' Apply lagtime to a regimen
-#'
-#' @param regimen PKPDsim regimen
-#' @param lagtime lagtime object, either single value / parameter name or vector of values/parameter names for all compartments.
-#' @param parameters parameter list, required if parameters are specified.
-#' @param cmt_mapping map of administration types to compartments, e.g. `list("oral" = 1, "infusion" = 2, "bolus" = 2)`.
-#'
-#' @export
-#' @return Original regimen with lagtime added to dose times
-apply_lagtime <- function(
-  regimen,
-  lagtime,
-  parameters,
-  cmt_mapping = NULL
-) {
-    if(class(lagtime) %in% c("numeric", "integer")) {
-      if(length(lagtime) == 1) {
-        regimen$dose_times <- regimen$dose_times + lagtime
-      } else {
-        regimen$dose_times <- regimen$dose_times + lagtime[regimen$cmt]
-      }
-    }
-    if(class(lagtime) %in% c("character")) {
-      if(length(lagtime) == 1) {
-        regimen$dose_times <- regimen$dose_times + parameters[[lagtime]]
-      } else {
-        if(is.null(regimen$cmt)) {
-          if(!is.null(cmt_mapping)) {
-            regimen$cmt <- as.numeric(cmt_mapping[regimen$type])
-          } else {
-            regimen$cmt <- rep(1, length(regimen$dose_times))
-          }
-        }
-        par_tmp <- parameters
-        par_tmp[["0"]] <- 0
-        regimen$dose_times <- regimen$dose_times + as.numeric(unlist(par_tmp[lagtime[regimen$cmt]]))
-      }
-    }
-    return(regimen)
-}
diff --git a/R/compile_sim_cpp.R b/R/compile_sim_cpp.R
index 417d75c1..32393448 100644
--- a/R/compile_sim_cpp.R
+++ b/R/compile_sim_cpp.R
@@ -270,10 +270,18 @@ compile_sim_cpp <- function(
   if(length(grep("-w", flg)) == 0) {
     Sys.setenv("PKG_CXXFLAGS" = paste(flg, "-w"))
   }
+
   if(compile) {
-    Rcpp::sourceCpp(code=sim_func, rebuild = TRUE, env = globalenv(), verbose = verbose, showOutput = verbose)
+    Rcpp::sourceCpp(
+      code = sim_func,
+      rebuild = TRUE, 
+      env = globalenv(), 
+      verbose = verbose, 
+      showOutput = verbose
+    )
     Sys.setenv("PKG_CXXFLAGS" = flg)
   }
+
   return(list(
     ode = ode_def_cpp,
     cpp = sim_func
diff --git a/R/create_event_table.R b/R/create_event_table.R
index 6567d0ba..4deb4fd5 100644
--- a/R/create_event_table.R
+++ b/R/create_event_table.R
@@ -152,10 +152,9 @@ create_event_table <- function(
   design <- dos[order(dos$t, -dos$dose),]
   if(!is.null(t_obs) && length(t_obs) != 0) { # make sure observation times are in dataset
     t_obs <- round(t_obs, 6)
-    t_diff <- setdiff(t_obs, design$t)
-    if(length(t_diff) > 0) {
-      design[(length(design[,1])+1) : (length(design[,1])+length(t_diff)),] <- cbind(
-         t = t_diff,
+    if(length(t_obs) > 0) {
+      design[(length(design[,1])+1) : (length(design[,1])+length(t_obs)),] <- cbind(
+         t = t_obs,
          dose = 0,
          type = 0,
          dum = 0,
diff --git a/R/parse_lagtime.R b/R/parse_lagtime.R
new file mode 100644
index 00000000..d87e88f2
--- /dev/null
+++ b/R/parse_lagtime.R
@@ -0,0 +1,27 @@
+#' Parse lagtime specified to main sim() function
+#' 
+#' @inheritParams sim
+#' 
+#' @returns a vector of character or numeric values
+#' 
+parse_lagtime <- function(
+  lagtime,
+  ode,
+  parameters
+) {
+  lagtime_ode <- attr(ode, "lagtime")
+  # override from ode if not specified by user and defined in ode
+  if(is.null(lagtime) && !is.null(lagtime_ode) && lagtime_ode[1] != "NULL" && lagtime_ode[1] != "undefined") {
+    lagtime <- lagtime_ode
+  }
+  if(is.null(lagtime)) {
+    lagtime <- 0
+  }
+  if(inherits(lagtime, "character")) {
+    idx <- grep("[a-zA-Z]", lagtime) # only pick character names, not "0"
+    if(! all(lagtime[idx] %in% names(parameters))) {
+      warning("Lagtime parameter(s) not found. Please check model and parameters.")
+    }
+  } 
+  lagtime
+}
\ No newline at end of file
diff --git a/R/sim.R b/R/sim.R
index f7dd0c50..8524f6a1 100644
--- a/R/sim.R
+++ b/R/sim.R
@@ -147,14 +147,7 @@ sim <- function (ode = NULL,
       regimen <- merge_regimen(list(regimen, regimen_dupl))
     }
   }
-  if(!is.null(attr(ode, "lagtime")) && attr(ode, "lagtime")[1] != "undefined" && attr(ode, "lagtime")[1] != "NULL") {
-    if(is.null(lagtime)) { # only override from metadata if not specified by user
-      lagtime <- attr(ode, "lagtime")
-    }
-  }
-  if(!is.null(lagtime)) {
-    regimen <- apply_lagtime(regimen, lagtime, parameters, attr(ode, "cmt_mapping"))
-  }
+  lagtime <- parse_lagtime(lagtime, ode, parameters)
   if(!is.null(attr(ode, "dose")$duration_scale)) {
     regimen <- apply_duration_scale(
       regimen,
@@ -234,6 +227,9 @@ sim <- function (ode = NULL,
     } else {
       size <- attr(analytical, "size")
     }
+    if(is.null(lagtime)) {
+      lagtime <- rep(0, size) # needs to have at least 1 zero value, cannot be NULL when passed to cpp func
+    }
     if(is.null(ode) && is.null(analytical)) {
       stop("Please specify at least the required arguments 'ode' or 'analytical' for simulations.")
     }
@@ -507,7 +503,7 @@ sim <- function (ode = NULL,
 
     #################### Main call to ODE solver / analytical eq solver #######################
     if(!is.null(ode)) {
-      tmp <- ode(A_init, design_i, p_i, iov_bins, int_step_size)
+      tmp <- ode(A_init, design_i, p_i, iov_bins, lagtime, int_step_size)
     } else {
       tmp <- analytical_eqn_wrapper(analytical, design_i, p_i)
     }
@@ -602,6 +598,12 @@ sim <- function (ode = NULL,
   all_names <- unique(c(par_names, cov_names, var_names))
   all_names <- intersect(all_names, names(comb)) # only cols that appear in data
 
+  ## remove dose-times from regimen
+  ## but leave ones that were also in t_obs
+  ## also leave extra obs when bolus, because this may be needed (controlled below using `extra_t_obs`)
+  dose_times <- design_i[design_i$evid > 0, ]$t
+  comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, comb$y, sep="_"))),]
+
   if(!extra_t_obs) {
     ## include the observations at which a bolus dose is added into the output object too
     comb <- comb[!duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, sep="_")),]
diff --git a/R/sim_core.R b/R/sim_core.R
index 44fe7bd8..6196940d 100644
--- a/R/sim_core.R
+++ b/R/sim_core.R
@@ -1,22 +1,28 @@
 #' Only core function of the simulation function, always just returns observations.
 #' Mostly useful for estimations / optimal design. Has no checks (for speed)!
 #'
+#' @inheritParams sim
 #' @param sim_object list with design and simulation parameters
-#' @param ode ode
-#' @param duplicate_t_obs allow duplicate t_obs in output? E.g. for optimal design calculations when t_obs = c(0,1,2,2,3). Default is FALSE.
-#' @param t_init time of initialization of the ODE system. Usually 0.
+#' 
 #' @export
-#' @return Data frame with simulation results
+#' 
+#' @return data.frame with simulation results
+#' 
 sim_core <- function(
   sim_object = NULL,
   ode,
   duplicate_t_obs = FALSE,
-  t_init = 0) {
-  tmp <- ode(A = sim_object$A_init,
-             design = sim_object$design,
-             par = sim_object$p,
-             iov_bins = sim_object$iov_bins,
-             step_size = sim_object$int_step_size)
+  t_init = 0,
+  lagtime = c(0)
+) {
+  tmp <- ode(
+    A = sim_object$A_init,
+    design = sim_object$design,
+    par = sim_object$p,
+    iov_bins = sim_object$iov_bins,
+    lagtime = lagtime,
+    step_size = sim_object$int_step_size
+  )
   out <- data.frame(t = tmp$time, y = tmp$obs, obs_type = tmp$obs_type)
   if(duplicate_t_obs) {
     # use match to ensure that duplicates in t_obs is possible
diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp
index ce8487a9..66ae6848 100644
--- a/inst/cpp/sim.cpp
+++ b/inst/cpp/sim.cpp
@@ -63,8 +63,97 @@ void pk_code (int i, std::vector<double> times, std::vector<double> doses, doubl
   // insert custom pk event code
 }
 
+NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) {
+   NumericVector lagtime_numeric;
+  if (TYPEOF(lagtime) == REALSXP) {
+    lagtime_numeric = as<NumericVector>(lagtime);
+  } else if (TYPEOF(lagtime) == STRSXP) {
+    CharacterVector lagtime_char = as<CharacterVector>(lagtime);
+    lagtime_numeric = NumericVector(lagtime_char.size());
+    for (int i = 0; i < lagtime_char.size(); i++) {
+      String param_name = lagtime_char[i];
+      if (parameters.containsElementNamed(param_name.get_cstring())) {
+        lagtime_numeric[i] = as<double>(parameters[param_name]);
+      } else {
+        lagtime_numeric[i] = 0.0; // default value
+      }
+    }
+  } else {
+    stop("lagtime must be either numeric or character vector");
+  }
+  return(lagtime_numeric);
+}
+
+List apply_lagtime(List design, NumericVector tlag, int n_comp) {
+
+  List new_design = clone(design);
+  std::vector<double> times = as<std::vector<double> >(new_design["t"]);
+  std::vector<int> evid = as<std::vector<int> >(new_design["evid"]);
+  std::vector<int> rate = as<std::vector<int> >(new_design["rate"]);
+  std::vector<int> cmt = as<std::vector<int> >(new_design["dose_cmt"]);
+  
+  NumericVector lagtime;
+  if(tlag.size() < n_comp) { // fill in with zeroes, if needed
+    int current_size = tlag.size();
+    lagtime = NumericVector(n_comp);
+    for(int i = 0; i < current_size; i++) {
+      lagtime[i] = tlag[i];
+    }
+    for(int i = current_size; i < n_comp; i++) {
+      lagtime[i] = 0.0;
+    }
+  } else {
+    lagtime = tlag;
+  }
+
+  // Apply lagtime to dose events (evid == 1)  
+  for(int i = 0; i < times.size(); i++) {
+    if(evid[i] == 1 | (evid[i] == 2 & rate[i] != 0)) { // dose events or infusion stop events
+      times[i] += lagtime[cmt[i]-1];
+    }
+  }
+
+  // Create sorted index according to "t"
+  std::vector<size_t> indices(times.size());
+  std::iota(indices.begin(), indices.end(), 0); // Fill with 0, 1, 2, ...
+  std::sort(indices.begin(), indices.end(), [&times](size_t i1, size_t i2) {
+    return times[i1] < times[i2];
+  });
+  
+  // Reorder all elements in `t` in new_design
+  std::vector<double> sorted_times(times.size());
+  for (size_t i = 0; i < indices.size(); i++) {
+    sorted_times[i] = times[indices[i]];
+  }
+  new_design["t"] = sorted_times;
+  
+  // Sort all other vectors in the design object
+  for (const char* key : {"dose", "type", "dum", "dose_cmt", "t_inf", "evid", "bioav", "rate", "obs_type"}) {
+    if (new_design.containsElementNamed(key)) {
+      SEXP vec = new_design[key];
+      if (TYPEOF(vec) == REALSXP) {
+        std::vector<double> old_vec = as<std::vector<double> >(vec);
+        std::vector<double> new_vec(old_vec.size());
+        for (size_t i = 0; i < indices.size(); i++) {
+          new_vec[i] = old_vec[indices[i]];
+        }
+        new_design[key] = new_vec;
+      } else if (TYPEOF(vec) == INTSXP) {
+        std::vector<int> old_vec = as<std::vector<int> >(vec);
+        std::vector<int> new_vec(old_vec.size());
+        for (size_t i = 0; i < indices.size(); i++) {
+          new_vec[i] = old_vec[indices[i]];
+        }
+        new_design[key] = new_vec;
+      }
+    }
+  }
+  
+  return new_design;
+}
+
 // [[Rcpp::export]]
-List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_bins, double step_size) {
+List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_bins, SEXP lagtime, double step_size) {
   std::vector<double> t;
   std::vector<state_type> y;
   // insert observation variable definition
@@ -72,14 +161,19 @@ List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_
   std::vector<double> times, doses, dummy, rates;
   std::vector<int> dose_cmt, dose_type, evid, obs_type, y_type;
   // insert variable definitions
-  times = as<std::vector<double> >(design["t"]);
-  doses = as<std::vector<double> >(design["dose"]);
-  evid = as<std::vector<int> >(design["evid"]);
-  dummy = as<std::vector<double> >(design["dum"]);
-  rates = as<std::vector<double> >(design["rate"]);
-  dose_cmt = as<std::vector<int> >(design["dose_cmt"]);
-  dose_type = as<std::vector<int> >(design["type"]);
-  obs_type = as<std::vector<int> >(design["obs_type"]);
+
+  // Handle lagtime parameter - can be numeric or character
+  NumericVector lagtime_numeric = lagtime_to_numeric(lagtime, par);
+  List events = apply_lagtime(design, lagtime_numeric, n_comp);
+  
+  times = as<std::vector<double> >(events["t"]);
+  doses = as<std::vector<double> >(events["dose"]);
+  evid = as<std::vector<int> >(events["evid"]);
+  dummy = as<std::vector<double> >(events["dum"]);
+  rates = as<std::vector<double> >(events["rate"]);
+  dose_cmt = as<std::vector<int> >(events["dose_cmt"]);
+  dose_type = as<std::vector<int> >(events["type"]);
+  obs_type = as<std::vector<int> >(events["obs_type"]);
   int len = times.size();
   int start;
   memset(rate, 0, sizeof(rate));
diff --git a/man/apply_lagtime.Rd b/man/apply_lagtime.Rd
deleted file mode 100644
index 3718dc5c..00000000
--- a/man/apply_lagtime.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/apply_lagtime.R
-\name{apply_lagtime}
-\alias{apply_lagtime}
-\title{Apply lagtime to a regimen}
-\usage{
-apply_lagtime(regimen, lagtime, parameters, cmt_mapping = NULL)
-}
-\arguments{
-\item{regimen}{PKPDsim regimen}
-
-\item{lagtime}{lagtime object, either single value / parameter name or vector of values/parameter names for all compartments.}
-
-\item{parameters}{parameter list, required if parameters are specified.}
-
-\item{cmt_mapping}{map of administration types to compartments, e.g. \code{list("oral" = 1, "infusion" = 2, "bolus" = 2)}.}
-}
-\value{
-Original regimen with lagtime added to dose times
-}
-\description{
-Apply lagtime to a regimen
-}
diff --git a/man/parse_lagtime.Rd b/man/parse_lagtime.Rd
new file mode 100644
index 00000000..658ec85d
--- /dev/null
+++ b/man/parse_lagtime.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/parse_lagtime.R
+\name{parse_lagtime}
+\alias{parse_lagtime}
+\title{Parse lagtime specified to main sim() function}
+\usage{
+parse_lagtime(lagtime, ode)
+}
+\arguments{
+\item{lagtime}{either a value (numeric) or a parameter (character) or NULL.}
+
+\item{ode}{function describing the ODE system}
+}
+\value{
+a vector of character or numeric values
+}
+\description{
+Parse lagtime specified to main sim() function
+}
diff --git a/man/sim_core.Rd b/man/sim_core.Rd
index 502bdb81..6e2d35e6 100644
--- a/man/sim_core.Rd
+++ b/man/sim_core.Rd
@@ -5,19 +5,27 @@
 \title{Only core function of the simulation function, always just returns observations.
 Mostly useful for estimations / optimal design. Has no checks (for speed)!}
 \usage{
-sim_core(sim_object = NULL, ode, duplicate_t_obs = FALSE, t_init = 0)
+sim_core(
+  sim_object = NULL,
+  ode,
+  duplicate_t_obs = FALSE,
+  t_init = 0,
+  lagtime = c(0)
+)
 }
 \arguments{
 \item{sim_object}{list with design and simulation parameters}
 
-\item{ode}{ode}
+\item{ode}{function describing the ODE system}
 
 \item{duplicate_t_obs}{allow duplicate t_obs in output? E.g. for optimal design calculations when t_obs = c(0,1,2,2,3). Default is FALSE.}
 
-\item{t_init}{time of initialization of the ODE system. Usually 0.}
+\item{t_init}{initialization time before first dose, default 0.}
+
+\item{lagtime}{either a value (numeric) or a parameter (character) or NULL.}
 }
 \value{
-Data frame with simulation results
+data.frame with simulation results
 }
 \description{
 Only core function of the simulation function, always just returns observations.
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 9dc6bb2c..906863f0 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -1,6 +1,16 @@
 mod_1cmt_iv <- new_ode_model("pk_1cmt_iv")
 mod_2cmt_iv <- new_ode_model("pk_2cmt_iv")
 mod_1cmt_oral <- new_ode_model("pk_1cmt_oral")
+mod_1cmt_oral_lagtime <- new_ode_model(
+  code = "
+    dAdt[0] = -KA * A[0]
+    dAdt[1] = +KA * A[0] -(CL/V) * A[1]
+  ",
+  lagtime = c("TLAG", 0),
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 1, bioav = 1),
+  parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
+)
 oral_1cmt_allometric <- new_ode_model( # also timevarying and dose-dependence factor
   code = "
       if(t<168.0) {
diff --git a/tests/testthat/test_apply_lagtime.R b/tests/testthat/test_apply_lagtime.R
deleted file mode 100644
index 5ff1c599..00000000
--- a/tests/testthat/test_apply_lagtime.R
+++ /dev/null
@@ -1,60 +0,0 @@
-reg1 <- new_regimen(
-  amt = 1000,
-  n = 12,
-  interval = 12,
-  type = "oral"
-)
-
-test_that("lagtime applied to regimens using parameters", {
-  lag1 <- apply_lagtime(
-    regimen = reg1,
-    lagtime = "TLAG",
-    parameters = list(CL = 5, V = 50, TLAG = .5)
-  )
-  expect_true("regimen" %in% class(lag1))
-  expect_equal(lag1$dose_times, reg1$dose_times + 0.5)
-})
-
-test_that("lagtime applied to regimens using lagtime arg", {
-  lag2 <- apply_lagtime(regimen = reg1, lagtime = .75)
-  expect_true("regimen" %in% class(lag2))
-  expect_equal(lag2$dose_times, reg1$dose_times + 0.75)
-})
-
-test_that("lagtime applied when multiple compartments and no compartment specified in regimen", {
-  lag3 <- apply_lagtime(
-    regimen = reg1,
-    lagtime = c("TLAG", 0, 0),
-    parameters = list(CL = 5, V = 50, TLAG = .5)
-  )
-  expect_true("regimen" %in% class(lag3))
-  expect_equal(lag3$dose_times, reg1$dose_times + 0.5)
-})
-
-test_that("lagtime applied when multiple compartments and no compartment specified in regimen but cmt_mapping available (oral)", {
-  lag4a <- apply_lagtime(
-    regimen = reg1,
-    lagtime = c("TLAG", 0, 0),
-    parameters = list(CL = 5, V = 50, TLAG = .5),
-    cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
-  )
-  expect_true("regimen" %in% class(lag4a))
-  expect_equal(lag4a$dose_times, reg1$dose_times + 0.5)
-})
-
-test_that("lagtime applied when multiple compartments and no compartment specified in regimen but cmt_mapping available (infusion): dose times should not be altered", {
-  reg2 <- new_regimen(
-    amt = 1000,
-    n = 12,
-    interval = 12,
-    type = "infusion"
-  )
-  lag4b <- apply_lagtime(
-    regimen = reg2,
-    lagtime = c("TLAG", 0, 0),
-    parameters = list(CL = 5, V = 50, TLAG = .5),
-    cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
-  )
-  expect_true("regimen" %in% class(lag4b))
-  expect_equal(lag4b$dose_times, reg2$dose_times)
-})
diff --git a/tests/testthat/test_is_positive_definite.R b/tests/testthat/test_is_positive_definite.R
index 4ccc82be..8dd334cd 100644
--- a/tests/testthat/test_is_positive_definite.R
+++ b/tests/testthat/test_is_positive_definite.R
@@ -105,7 +105,9 @@ test_that("Matrix where numerical solver issues result in very tiny complex comp
   ),
   nrow = 16
   )
-  is_x_positive_def <- is_positive_definite(x)
+  suppressWarnings(
+    is_x_positive_def <- is_positive_definite(x)
+  )
   expect_true(is_x_positive_def)
 })
 
diff --git a/tests/testthat/test_parse_lagtime.R b/tests/testthat/test_parse_lagtime.R
new file mode 100644
index 00000000..43dfb03b
--- /dev/null
+++ b/tests/testthat/test_parse_lagtime.R
@@ -0,0 +1,83 @@
+test_that("parse_lagtime handles NULL lagtime input", {
+  # Mock ODE object without lagtime attribute
+  ode <- list()
+  parameters <- list(CL = 5, V = 50)
+  
+  result <- parse_lagtime(NULL, ode, parameters)
+  expect_equal(result, 0)
+})
+
+test_that("parse_lagtime uses ODE lagtime when user lagtime is NULL", {
+  # Mock ODE object with lagtime attribute
+  ode <- list()
+  attr(ode, "lagtime") <- "TLAG"
+  parameters <- list(CL = 5, V = 50, TLAG = 0.6)
+  
+  result <- parse_lagtime(NULL, ode, parameters)
+  expect_equal(result, "TLAG")
+})
+
+test_that("parse_lagtime ignores ODE lagtime when user provides lagtime", {
+  # Mock ODE object with lagtime attribute
+  ode <- list()
+  attr(ode, "lagtime") <- "TLAG"
+  parameters <- list(CL = 5, V = 50, TLAG = 0.6)
+  
+  result <- parse_lagtime(1.5, ode, parameters)
+  expect_equal(result, 1.5)
+})
+
+test_that("parse_lagtime handles character lagtime with valid parameters", {
+  ode <- list()
+  parameters <- list(CL = 5, V = 50, TLAG = 0.6)
+  lagtime <- c("TLAG", "0", "0")
+  
+  result <- parse_lagtime(lagtime, ode, parameters)
+  expect_equal(result, c("TLAG", "0", "0"))
+})
+
+test_that("parse_lagtime warns when character lagtime parameters not found", {
+  ode <- list()
+  parameters <- list(CL = 5, V = 50)
+  lagtime <- c("TLAG", "0", "0")
+  
+  expect_warning(
+    parse_lagtime(lagtime, ode, parameters),
+    "Lagtime parameter\\(s\\) not found. Please check model and parameters."
+  )
+})
+
+test_that("parse_lagtime handles numeric lagtime", {
+  ode <- list()
+  parameters <- list(CL = 5, V = 50)
+  lagtime <- 1.5
+  
+  result <- parse_lagtime(lagtime, ode, parameters)
+  expect_equal(result, 1.5)
+})
+
+test_that("parse_lagtime ignores NULL and undefined ODE lagtime", {
+  # Test NULL lagtime in ODE
+  ode1 <- list()
+  attr(ode1, "lagtime") <- "NULL"
+  parameters <- list(CL = 5, V = 50)
+  
+  result1 <- parse_lagtime(NULL, ode1, parameters)
+  expect_equal(result1, 0)
+  
+  # Test undefined lagtime in ODE
+  ode2 <- list()
+  attr(ode2, "lagtime") <- "undefined"
+  
+  result2 <- parse_lagtime(NULL, ode2, parameters)
+  expect_equal(result2, 0)
+})
+
+test_that("parse_lagtime handles mixed character and numeric lagtime", {
+  ode <- list()
+  parameters <- list(CL = 5, V = 50, TLAG1 = 0.6, TLAG2 = 1.2)
+  lagtime <- c("TLAG1", "0", "TLAG2")
+  
+  result <- parse_lagtime(lagtime, ode, parameters)
+  expect_equal(result, c("TLAG1", "0", "TLAG2"))
+})
diff --git a/tests/testthat/test_sim.R b/tests/testthat/test_sim.R
index 0e3892b7..4b488bec 100644
--- a/tests/testthat/test_sim.R
+++ b/tests/testthat/test_sim.R
@@ -20,19 +20,21 @@ test_that("return_event_table=TRUE returns an appropriate event table", {
   )
   expect_equal(
     evtab1,
-    data.frame(
-      t = c(0, 2, 12, 14, 24, 26, 48, 48),
-      dose = c(100, 0, 100, 0, 100, 0, 0, 0),
-      type = c(1, 1, 1, 1, 1, 1, 0, 0),
-      dum = c(0, 1, 0, 1, 0, 1, 0, 0),
-      dose_cmt = c(1, 1, 1, 1, 1, 1, 0, 0),
-      t_inf = c(2, 0, 2, 0, 2, 0, 0, 0),
-      evid = c(1, 2, 1, 2, 1, 2, 0, 0),
-      bioav = c(1, 0, 1, 0, 1, 0, 0, 0),
-      rate = c(50,-50, 50,-50, 50,-50, 0, 0),
-      obs_type = c(0, 1, 0, 0, 0, 0, 1, 1)
-    )
-  )
+    structure(list(
+      t = c(0, 2, 2, 12, 14, 24, 26, 48, 48), 
+      dose = c(100, 0, 0, 100, 0, 100, 0, 0, 0), 
+      type = c(1, 0, 1, 1, 1, 1, 1, 0, 0), 
+      dum = c(0, 0, 1, 0, 1, 0, 1, 0, 0), 
+      dose_cmt = c(1, 0, 1, 1, 1, 1, 1, 0, 0), 
+      t_inf = c(2, 0, 0, 2, 0, 2, 0, 0, 0), 
+      evid = c(1, 0, 2, 1, 2, 1, 2, 0, 0), 
+      bioav = c(1, 0, 0, 1, 0, 1, 0, 0, 0), 
+      rate = c(50, 0, -50, 50, -50, 50, -50, 0, 0), 
+      obs_type = c(0, 1, 1, 0, 0, 0, 0, 1, 1)
+    ), 
+    row.names = c(1L, 3L, 2L, 4L, 5L, 6L, 7L, 8L, 9L), 
+    class = "data.frame"
+  ))
 })
 
 test_that("return_event_table=TRUE returns an appropriate event table with covariate", {
@@ -47,30 +49,27 @@ test_that("return_event_table=TRUE returns an appropriate event table with covar
   )
   expect_equal(
     evtab2,
-    structure(
-      list(
-        t = c(0, 0, 2, 12, 14, 24, 24, 26, 48),
-        dose = c(0,  100, 0, 100, 0, 0, 100, 0, 0),
-        type = c(0, 1, 1, 1, 1, 0, 1, 1, 0),
-        dum = c(0, 0, 1, 0, 1, 0, 0, 1, 0),
-        dose_cmt = c(0, 1, 1, 1, 1, 0, 1, 1, 0),
-        t_inf = c(0, 2, 0, 2, 0, 0, 2, 0, 0),
-        evid = c(2, 1, 2, 1, 2, 2, 1, 2, 0),
-        bioav = c(1, 1, 0, 1, 0, 1, 1, 0, 0),
-        rate = c(0, 50,-50, 50,-50, 0, 50,-50, 0),
-        cov_CRCL = c(70, 70, 70, 70, 70, 80, 80, 80, 80),
-        cov_t_CRCL = c(0, 0, 0, 0, 0, 24, 24, 24, 24),
-        gradients_CRCL = c(0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0, 0, 0, 0),
-        cov_WT = c(70, 70, 70, 70, 70, 70, 70, 70, 70),
-        cov_t_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0),
-        gradients_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0),
-        obs_type = c(0, 0, 1, 0, 0, 0, 0, 0, 1)
-      ),
-      row.names = c(2L, 1L, 3L, 4L, 5L,
-                    6L, 7L, 8L, 9L),
-      class = "data.frame"
-    )
-  )
+    structure(list(
+      t = c(0, 0, 2, 2, 12, 14, 24, 24, 26, 48), 
+      dose = c(0, 100, 0, 0, 100, 0, 0, 100, 0, 0), 
+      type = c(0, 1, 0, 1, 1, 1, 0, 1, 1, 0), 
+      dum = c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0), 
+      dose_cmt = c(0, 1, 0, 1, 1, 1, 0, 1, 1, 0), 
+      t_inf = c(0, 2, 0, 0, 2, 0, 0, 2, 0, 0), 
+      evid = c(2, 1, 0, 2, 1, 2, 2, 1, 2, 0), 
+      bioav = c(1, 1, 0, 0, 1, 0, 1, 1, 0, 0), 
+      rate = c(0, 50, 0, -50, 50, -50, 0, 50, -50, 0), 
+      cov_CRCL = c(70, 70, 70, 70, 70, 70, 80, 80, 80, 80), 
+      cov_t_CRCL = c(0, 0, 0, 0, 0, 0, 24, 24, 24, 24), 
+      gradients_CRCL = c(0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0, 0, 0, 0), 
+      cov_WT = c(70, 70, 70, 70, 70, 70, 70, 70, 70, 70), 
+      cov_t_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 
+      gradients_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 
+      obs_type = c(0, 0, 1, 1, 0, 0, 0, 0, 0, 1)
+    ), 
+    row.names = c(2L, 1L, 4L, 3L, 5L, 6L, 8L, 7L, 9L, 10L), 
+    class = "data.frame"
+  ))
 })
 
 test_that("sim works properly for a model where bioavailability is dependent on dose", {
@@ -183,7 +182,6 @@ test_that("covariates and doses are shifted correctly when t_init != 0", {
 })
 
 
-
 test_that("covariates_table and doses are shifted correctly when t_init != 0", {
   cov_table <- data.frame(
     id = c(1, 1),
diff --git a/tests/testthat/test_sim_core.R b/tests/testthat/test_sim_core.R
index ce09b429..fd81ca48 100644
--- a/tests/testthat/test_sim_core.R
+++ b/tests/testthat/test_sim_core.R
@@ -16,3 +16,36 @@ test_that("sim core works", {
   }
   expect_equal(f1(), f2())
 })
+
+test_that("sim core works for absorption model with lagtime", {  
+  reg <- new_regimen(amt = 100, n = 5, interval = 12, t_inf = 1, type = "infusion")
+  par <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
+
+  ## have to be explicit about t_obs with sim_core!
+  f_ref <- function() {
+    res <- sim(
+      ode = mod_1cmt_oral_lagtime, 
+      regimen = reg, 
+      parameters = par, 
+      only_obs = TRUE, 
+      t_obs=c(0:24)
+    )$y
+    return(res)
+  }
+  f_core <- function() {
+    obj <- sim(
+      ode = mod_1cmt_oral_lagtime, 
+      regimen = reg, 
+      parameters = par, 
+      only_obs = TRUE, 
+      t_obs=c(0:24), 
+      return_design=TRUE
+    )
+    sim_core(
+      obj, 
+      ode = mod_1cmt_oral_lagtime, 
+      lagtime = c(0.83, 0)
+    )$y ## Lagtime parameter needed!
+  }
+  expect_equal(f_ref(), f_core())
+})
diff --git a/tests/testthat/test_sim_lagtime.R b/tests/testthat/test_sim_lagtime.R
index 91951509..77a8c5f7 100644
--- a/tests/testthat/test_sim_lagtime.R
+++ b/tests/testthat/test_sim_lagtime.R
@@ -2,21 +2,25 @@ test_that("dose dump after lagtime in correct order in output data", {
   skip_on_cran()
   reg <- new_regimen(amt = 500, n = 4, interval = 12, type = 'oral')
   pars <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
-  mod <- new_ode_model(
-    code = "
-      dAdt[0] = -KA * A[0]
-      dAdt[1] = +KA * A[0] -(CL/V) * A[1]
-  ",
-    lagtime = c("TLAG", 0),
-    obs = list(cmt = 2, scale = "V"),
-    dose = list(cmt = 1, bioav = 1),
-    parameters = pars
-  )
   dat <- sim_ode(
-    ode = mod,
+    ode = mod_1cmt_oral_lagtime,
     regimen = reg,
     parameters = pars,
     only_obs = FALSE
   )
-  expect_equal(round(dat[dat$t == 12.83 & dat$comp == 1,]$y, 1), c(1.2, 501.2))
+  ## Change after RXR-2394: time of TLAG not in dataset anymore unless requested by user in t_obs
+  ## Before change: expect_equal(round(dat[dat$t == 12.83 & dat$comp == 1,]$y, 1), c(1.2, 501.2))
+  ## After change:
+  expect_equal(nrow(dat[dat$t == 12.83,]), 0)
+  ## When grid requested by user, lagtime should be visible
+  dat <- sim_ode(
+    ode = mod_1cmt_oral_lagtime,
+    regimen = reg,
+    parameters = pars,
+    t_obs = seq(0, 1, 0.01),
+    only_obs = TRUE
+  )
+  tmp <- dat[dat$t >= 0.82 & dat$t <= 0.85, ]
+  expect_equal(tmp$t, c(0.82, 0.83, 0.84, 0.85))
+  expect_equal(round(tmp$y, 2), c(0, 0, 0.05, 0.10))
 })