diff --git a/R/RcppExports.R b/R/RcppExports.R
index f648dd00..87cb2eab 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -25,8 +25,8 @@ C_cohort_ctstm_sim <- function(R_CtstmTrans, R_CostsStateVals, R_QALYsStateVal,
     .Call('_hesim_C_cohort_ctstm_sim', PACKAGE = 'hesim', R_CtstmTrans, R_CostsStateVals, R_QALYsStateVal, live_states, start_state, start_age, times, clock, transition_types, progress, dr_qalys, dr_costs, type, zero_tol, abs_tol, rel_tol)
 }
 
-runTestODE <- function(p0, times, discount_rate = 0.0) {
-    .Call('_hesim_runTestODE', PACKAGE = 'hesim', p0, times, discount_rate)
+runCohortCtstmTestODE <- function(start_state, times, discount_rate = 0.0) {
+    .Call('_hesim_runCohortCtstmTestODE', PACKAGE = 'hesim', start_state, times, discount_rate)
 }
 
 tmax_max <- function(m) {
diff --git a/R/ctstm.R b/R/ctstm.R
index 2d88cd12..9495c0d1 100644
--- a/R/ctstm.R
+++ b/R/ctstm.R
@@ -808,7 +808,11 @@ IndivCtstm <- R6::R6Class("IndivCtstm",
 #' clock-reset (essentially fixed times in state), whereas in
 #' `CohortCtstm` the point masses are always clock-forward (Markov).
 #' 
-#' There are currently no relevant vignettes. 
+#' There is currently one relevant vignette.
+#' `vignette("markov-inhomogeneous-cohort-continuous")` shows how the
+#' individual-based time inhomogeneous Markov model in
+#' `vignette("markov-inhomogeneous-indiv")` can be developed as a
+#' cohort-based simulation.
 #' @name CohortCtstm
 NULL
 
@@ -965,3 +969,30 @@ CohortCtstm <- R6::R6Class("CohortCtstm",
     }
   ) # end public
 ) # end class
+
+#' @description 
+#' @internal
+run_CohortCtstmTestODE <- function(start_state=1, times=c(0,1,2,10), discount_rate=0.03) {
+    self <- runCohortCtstmTestODE(start_state - 1L, times, discount_rate)
+    ## NB: grp_id is already 1-based -- so do not update
+    self$stateprobs_ <-
+        as.data.table(self$stateprobs_)[,`:=`(sample=sample+1L,strategy_id=strategy_id+1L,
+                                             patient_id=patient_id+1L,state_id=state_id+1L)]
+    ev <- as.data.table(self$ev_)
+    self$qalys_ <- ev[outcome=="qaly"][
+       ,`:=`(sample=sample+1L,strategy_id=strategy_id+1L,patient_id=patient_id+1L,
+             state_id=state_id+1,category='qalys',qalys=value,value=NULL,outcome=NULL)]
+    self$qalys_$ly = ev[outcome=="ly","value"]
+    self$costs_ <- do.call(rbind,
+                           lapply(unique(grep("Category",ev$outcome,value=TRUE)),
+                                  function(pattern) {
+                                      ev2 <- ev[outcome==pattern][
+                                         ,`:=`(sample=sample+1L,strategy_id=strategy_id+1L,
+                                               patient_id=patient_id+1L,state_id=state_id+1L,
+                                               category=pattern,costs=value,value=NULL,
+                                               outcome=NULL)]
+                                      ev2
+                                  }))
+    self$ev_ = NULL
+    self
+}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index fbf409f8..9bf0dea0 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -115,16 +115,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// runTestODE
-Rcpp::List runTestODE(arma::vec p0, arma::vec times, double discount_rate);
-RcppExport SEXP _hesim_runTestODE(SEXP p0SEXP, SEXP timesSEXP, SEXP discount_rateSEXP) {
+// runCohortCtstmTestODE
+Rcpp::List runCohortCtstmTestODE(int start_state, arma::vec times, double discount_rate);
+RcppExport SEXP _hesim_runCohortCtstmTestODE(SEXP start_stateSEXP, SEXP timesSEXP, SEXP discount_rateSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< arma::vec >::type p0(p0SEXP);
+    Rcpp::traits::input_parameter< int >::type start_state(start_stateSEXP);
     Rcpp::traits::input_parameter< arma::vec >::type times(timesSEXP);
     Rcpp::traits::input_parameter< double >::type discount_rate(discount_rateSEXP);
-    rcpp_result_gen = Rcpp::wrap(runTestODE(p0, times, discount_rate));
+    rcpp_result_gen = Rcpp::wrap(runCohortCtstmTestODE(start_state, times, discount_rate));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -706,7 +706,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_hesim_C_mce", (DL_FUNC) &_hesim_C_mce, 6},
     {"_hesim_C_enmbpi", (DL_FUNC) &_hesim_C_enmbpi, 6},
     {"_hesim_C_cohort_ctstm_sim", (DL_FUNC) &_hesim_C_cohort_ctstm_sim, 16},
-    {"_hesim_runTestODE", (DL_FUNC) &_hesim_runTestODE, 3},
+    {"_hesim_runCohortCtstmTestODE", (DL_FUNC) &_hesim_runCohortCtstmTestODE, 3},
     {"_hesim_tmax_max", (DL_FUNC) &_hesim_tmax_max, 1},
     {"_hesim_C_ctstm_summary", (DL_FUNC) &_hesim_C_ctstm_summary, 3},
     {"_hesim_C_rgengamma", (DL_FUNC) &_hesim_C_rgengamma, 4},
diff --git a/src/cohort-ctstm.cpp b/src/cohort-ctstm.cpp
index 731b5cf4..e77f6b91 100644
--- a/src/cohort-ctstm.cpp
+++ b/src/cohort-ctstm.cpp
@@ -36,14 +36,6 @@ namespace boost {
   }
 } // namespace boost::numeric::odeint
 
-#define CTSTM_OUT(OBJECT,COUNTER)	     \
-  OBJECT.sample_[COUNTER] = s;		     \
-  OBJECT.strategy_id_[COUNTER] = k;	     \
-  OBJECT.patient_id_[COUNTER] = i;					\
-  OBJECT.grp_id_[COUNTER] = transmod->obs_index_.get_grp_id();		\
-  OBJECT.patient_wt_[COUNTER] = transmod->obs_index_.get_patient_wt();	\
-  OBJECT.state_id_[COUNTER] = h;
-
 /***************************************************************************//** 
  * @ingroup ctstm
  * Simulate disease progression (i.e., a path through a multi-state model),
@@ -315,20 +307,33 @@ Rcpp::List C_cohort_ctstm_sim(Rcpp::Environment R_CtstmTrans,
 	  } // end j full_times_ loop
 	  for (int h = 0; h < n_states; ++h){
 	    for (int ti = 0; ti < n_times; ++ti){
-	      CTSTM_OUT(out,counter);
+	      out.sample_[counter] = s;
+	      out.strategy_id_[counter] = k;
+	      out.patient_id_[counter] = i;
+	      out.grp_id_[counter] = transmod->obs_index_.get_grp_id();
+	      out.patient_wt_[counter] = transmod->obs_index_.get_patient_wt();
+	      out.state_id_[counter] = h;
 	      out.t_[counter] = times[ti];
 	      out.prob_[counter] = report(ti, h);
 	      ++counter;                   
 	    } // end cycles loop
 	    // report for life-years
-	    CTSTM_OUT(out2,counter2);
+	    auto update_ev_out = [&](hesim::ev_out & object,int counter) {
+	      object.sample_[counter] = s;
+	      object.strategy_id_[counter] = k;
+	      object.patient_id_[counter] = i;
+	      object.grp_id_[counter] = transmod->obs_index_.get_grp_id();
+	      object.patient_wt_[counter] = transmod->obs_index_.get_patient_wt();
+	      object.state_id_[counter] = h;
+	    };
+	    update_ev_out(out2,counter2);
 	    out2.dr_[counter2] = 0.0; // assume no discounting for life-years
 	    out2.outcome_[counter2] = "ly";
 	    out2.value_[counter2] = report(n_times-1,n_states+h);
 	    ++counter2;
 	    // report for qalys
 	    if (valid_R_QALYsStateVal) {
-	      CTSTM_OUT(out2,counter2);
+	      update_ev_out(out2,counter2);
 	      out2.dr_[counter2] = dr_qalys;
 	      out2.outcome_[counter2] = "qaly";
 	      out2.value_[counter2] = report(n_times-1,2*n_states+h);
@@ -336,7 +341,7 @@ Rcpp::List C_cohort_ctstm_sim(Rcpp::Environment R_CtstmTrans,
 	    }
 	    // report for costs
 	    for (int cost_=0; cost_<n_costs; ++cost_) {
-	      CTSTM_OUT(out2,counter2);
+	      update_ev_out(out2,counter2);
 	      out2.dr_[counter2] = dr_costs;
 	      out2.outcome_[counter2] = "Category " + std::to_string(cost_+1);
 	      out2.value_[counter2] = report(n_times-1,(3+cost_)*n_states+h);
@@ -352,18 +357,17 @@ Rcpp::List C_cohort_ctstm_sim(Rcpp::Environment R_CtstmTrans,
 		      _["ev"]=out2.create_R_data_frame()));
 }
 
-#undef CTSTM_OUT
-
 // test example that does not use hesim
-class TestODE {
+class CohortCtstmTestODE {
 public:
   typedef arma::vec state_type;
   size_t n_states;
   double discount_rate;
-  TestODE(double discount_rate = 0.0)
+  CohortCtstmTestODE(double discount_rate = 0.0)
     : n_states(4), discount_rate(discount_rate) { }
   inline double hweibull(const double t, const double shape, const double scale) {
-    return R::dweibull(std::max(1e-100,t),shape,scale,0)/R::pweibull(t,shape,scale,0,0);
+    return R::dweibull(std::max(1e-100,t),shape,scale,0)/
+      R::pweibull(std::max(1e-100,t),shape,scale,0,0);
   }
   void operator() (const state_type &Y , state_type &dYdt, const double t)
   {
@@ -402,13 +406,17 @@ class TestODE {
     for (int j=0; j<4; ++j)
       dYdt(to[j] + 5*n_states) += starting_costs[to[j]] * Y(from[j]) * rates(j) * dr;
   }
-  arma::mat run(arma::vec p0, arma::vec times) {
+  Rcpp::List run(int start_state = 0, arma::vec times = {0.0, 1.0, 2.0, 10.0},
+		 double discount_rate = 0.0) {
+    this->discount_rate = discount_rate;
     using namespace boost::numeric::odeint;
     size_t n_times = times.size();
+    arma::vec p0(n_states*6); // probs, lys, qalys, 3*costs
+    p0(start_state) = 1.0;
     // combine the results
     arma::mat combined(n_times,p0.size());
     combined.row(0) = p0.t();
-    auto stepper = make_dense_output(1.0e-10, 1.0e-10, runge_kutta_dopri5<state_type>());
+    auto stepper = make_dense_output(1.0e-6, 1.0e-6, runge_kutta_dopri5<state_type>());
     for (size_t step_=1; step_<n_times; step_++) {
       size_t n = integrate_adaptive(stepper,
 				    [this](const state_type &Y , state_type &dYdt, const double t) {
@@ -420,12 +428,62 @@ class TestODE {
 				    times[step_]-times[step_-1]);
       combined.row(step_) = p0.t();
     }
-    return combined;
+    int counter=0;
+    hesim::stateprobs_out out(n_times * n_states);
+    for (int h = 0; h < n_states; ++h){
+      for (int ti = 0; ti < n_times; ++ti){
+	out.sample_[counter] = 0; // C-based indexing
+	out.strategy_id_[counter] = 0;
+	out.patient_id_[counter] = 0;
+	out.grp_id_[counter] = 1; // ???
+	out.patient_wt_[counter] = 1.0;
+	out.state_id_[counter] = h;
+	out.t_[counter] = times[ti];
+	out.prob_[counter] = combined(ti, h);
+	++counter;                   
+      } // end cycles loop
+    } // end 
+    hesim::ev_out out2(n_states * 5);
+    int counter2 = 0;
+    for (int h = 0; h < n_states; ++h) {
+      auto update_ev_out = [&](hesim::ev_out & object,int counter) {
+	object.sample_[counter] = 0;
+	object.strategy_id_[counter] = 0;
+	object.patient_id_[counter] = 0;
+	object.grp_id_[counter] = 1;
+	object.patient_wt_[counter] = 1.0;
+	object.state_id_[counter] = h;
+	object.dr_[counter] = discount_rate;
+      };
+      update_ev_out(out2,counter2);
+      out2.dr_[counter2] = 0.0;
+      out2.outcome_[counter2] = "ly";
+      out2.value_[counter2] = combined(n_times-1,1*n_states+h);
+      ++counter2;
+      update_ev_out(out2,counter2);
+      out2.outcome_[counter2] = "qaly";
+      out2.value_[counter2] = combined(n_times-1,2*n_states+h);
+      ++counter2;
+      update_ev_out(out2,counter2);
+      out2.outcome_[counter2] = "Category 1";
+      out2.value_[counter2] = combined(n_times-1,3*n_states+h);
+      ++counter2;
+      update_ev_out(out2,counter2);
+      out2.outcome_[counter2] = "Category 2";
+      out2.value_[counter2] = combined(n_times-1,4*n_states+h);
+      ++counter2;
+      update_ev_out(out2,counter2);
+      out2.outcome_[counter2] = "Category 3";
+      out2.value_[counter2] = combined(n_times-1,5*n_states+h);
+      ++counter2;
+    }
+    using namespace Rcpp;
+    return List::create(_["stateprobs_"]=out.create_R_data_frame(),
+			_["ev_"]=out2.create_R_data_frame());
   }
 };
 // [[Rcpp::export]]
-Rcpp::List runTestODE(arma::vec p0, arma::vec times, double discount_rate = 0.0) {
-  using namespace Rcpp;
-  return List::create(_("times")=times,
-		      _("Y")=TestODE(discount_rate).run(p0,times));
+Rcpp::List runCohortCtstmTestODE(int start_state, arma::vec times,
+				 double discount_rate = 0.0) {
+  return CohortCtstmTestODE().run(start_state,times,discount_rate);
 }
diff --git a/tests/testthat/test-ctstm.R b/tests/testthat/test-ctstm.R
index 280a6844..4242549d 100644
--- a/tests/testthat/test-ctstm.R
+++ b/tests/testthat/test-ctstm.R
@@ -595,6 +595,104 @@ test_that("Simulate costs and QALYs", {
   expect_true("ceac" %in% names(cea_pw))
 })
 
+test_that("Simulate a cohort-based continuous-time Markov model", {
+    strategies <- data.frame(strategy_id = 1L)
+    patients <- data.frame(patient_id = 1:2, intercept=1)
+    ## Multi-state model with transition specific models (as per testODE in the C++ code)
+    tmat <- rbind(c(NA, 1, NA, NA),
+                  c(2, NA, 3, NA),
+                  c(NA, NA, NA, 4),
+                  c(NA, NA, NA, NA))
+    ## short utility function
+    make_param <- function(shape,scale,dist="weibull")
+        params_surv(
+            coefs = list(
+                shape = data.frame(
+                    intercept = log(shape)
+                ),
+                scale = data.frame(
+                    intercept = log(scale)
+                )
+            ),
+            dist = dist)
+    fits = params_surv_list(make_param(0.8, 1.0), make_param(0.8, 1.0),
+                            make_param(0.8, 2.0), make_param(0.8, 3.0))
+    ## Simulation model
+    hesim_dat <- hesim_data(strategies = strategies,
+                            patients = patients)
+    fits_data <- expand(hesim_dat)
+    ## create_IndivCtstmTrans can also be used for cCTSTM:)
+    transmod <- create_IndivCtstmTrans(fits, input_data = fits_data,
+                                       trans_mat = tmat,
+                                       start_age=50,
+                                       start_state=1,
+                                       clock="mixt",
+                                       transition_types=c("time","time","time","time"))
+    ## utilities
+    utilities_tbl = stateval_tbl(data.frame(state_id=1:3,
+                                            est=c(0.9, 0.8, 0.7)),
+                                 dist="fixed")
+    utilities = create_StateVals(utilities_tbl, hesim_data=hesim_dat, n=1)
+    ## costs
+    costs_wlos_tbl = stateval_tbl(data.frame(state_id=1:3,
+                                             est=c(10000.0, 20000.0, 30000.0)),
+                                  dist="fixed")
+    costs_wlos = create_StateVals(costs_wlos_tbl, hesim_data=hesim_dat, n=1)
+    ## for the transition costs, "state_id" is actually the id for the transition
+    costs_transition_tbl <- stateval_tbl(tbl = data.frame(state_id = 1:4,
+                                                          est = c(10, 15, 20, 30)),
+                                         dist = "fixed")
+    costs_transition <- create_StateVals(costs_transition_tbl, n = 1,
+                                         time_reset = FALSE, method = "transition",
+                                         hesim_data = hesim_dat) 
+    costs_starting_tbl = stateval_tbl(data.frame(state_id=1:3,
+                                                 est=c(0.0, 2000.0, 3000.0)),
+                                      dist="fixed")
+    costs_starting = create_StateVals(costs_starting_tbl, hesim_data=hesim_dat, n=1,
+                                      method="starting")
+
+    test <- CohortCtstm$new(trans_model = transmod,
+                        utility_model = utilities,
+                        cost_models = list(costs_wlos,costs_transition,costs_starting))
+
+    ## 0.03 discount rate
+    test$sim(t=c(0,1,2,10), dr_qalys=0.03, dr_costs=0.03)
+    test2 <- hesim:::run_CohortCtstmTestODE(1, c(0,1,2,10), 0.03)
+    expect_true(max(abs(subset(test$stateprobs_, patient_id==1) - test2$stateprobs_)) < 5e-15)
+    expect_true(all(subset(test$qalys_, patient_id==1)[,-(9:10)] == test2$qalys_[,-(9:10)]))
+    expect_true(max(abs(subset(test$qalys_, patient_id==1)[,-8] - test2$qalys_[,-8])) < 5e-15)
+    expect_true(all(subset(test$costs_, patient_id==1)[,-9] == test2$costs_[,-9]))
+    expect_true(max(abs(subset(test$qalys_, patient_id==1)[,-8] - test2$qalys_[,-8])) < 5e-15)
+
+    ## zero discount
+    test$sim(t=c(0,1,2,10), dr_qalys=0.0, dr_costs=0.0)
+    test2 <- hesim:::run_CohortCtstmTestODE(1, c(0,1,2,10), 0.0)
+    expect_true(max(abs(subset(test$stateprobs_, patient_id==1) - test2$stateprobs_)) < 5e-15)
+    expect_true(all(subset(test$qalys_, patient_id==1)[,-(9:10)] == test2$qalys_[,-(9:10)]))
+    expect_true(max(abs(subset(test$qalys_, patient_id==1)[,-8] - test2$qalys_[,-8])) < 5e-15)
+    expect_true(all(subset(test$costs_, patient_id==1)[,-9] == test2$costs_[,-9]))
+    expect_true(max(abs(subset(test$qalys_, patient_id==1)[,-8] - test2$qalys_[,-8])) < 5e-15)
+
+    ## start_state=2
+    transmod <- create_IndivCtstmTrans(fits, input_data = fits_data,
+                                   trans_mat = tmat,
+                                   start_age=50,
+                                   start_state=2,
+                                   clock="mixt",
+                                   transition_types=c("time","time","time","time"))
+    test <- CohortCtstm$new(trans_model = transmod,
+                            utility_model = utilities,
+                            cost_models = list(costs_wlos,costs_transition,costs_starting))
+    test$sim(t=c(0,1,2,10), dr_qalys=0.03, dr_costs=0.03)
+    test2 <- hesim:::run_CohortCtstmTestODE(2, c(0,1,2,10), 0.03)
+    expect_true(max(abs(subset(test$stateprobs_, patient_id==1) - test2$stateprobs_)) < 5e-15)
+    expect_true(all(subset(test$qalys_, patient_id==1)[,-(9:10)] == test2$qalys_[,-(9:10)]))
+    expect_true(max(abs(subset(test$qalys_, patient_id==1)[,-8] - test2$qalys_[,-8])) < 5e-15)
+    expect_true(all(subset(test$costs_, patient_id==1)[,-9] == test2$costs_[,-9]))
+    expect_true(max(abs(subset(test$qalys_, patient_id==1)[,-8] - test2$qalys_[,-8])) < 5e-15)
+})
+
+
 ## With a joint survival model
 test_that("IndivCtstm - joint", {
   # Clock-reset