diff --git a/NEWS.md b/NEWS.md
index 8fa3f358..9d0d5460 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,14 @@
 * Fix bug where if/else clauses in the model could cause the model to
   not predict the values correctly.
 
+* Fix bug so that `shrinkage()` calculation works
+
+* Fix bug so that you can mix 2 different `PopED` data bases in an
+  analysis without crashing R.  While this didn't occur with every
+  database clash, it more frequently occurred when you interleaved
+  `PopED` code between two different `PopED` databases, like in issue
+  #131.
+
 # babelmixr2 0.1.4
 
 * Added experimental `PopED` integration
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 97659764..29584c78 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -145,8 +145,8 @@ popedFree <- function() {
     .Call(`_babelmixr2_popedFree`)
 }
 
-popedSetup <- function(e, full) {
-    .Call(`_babelmixr2_popedSetup`, e, full)
+popedSetup <- function(e, eglobal, full) {
+    .Call(`_babelmixr2_popedSetup`, e, eglobal, full)
 }
 
 popedSolveIdME <- function(theta, umt, mt, ms, nend, id, totn) {
diff --git a/R/poped.R b/R/poped.R
index dcc8f317..81d36423 100644
--- a/R/poped.R
+++ b/R/poped.R
@@ -17,13 +17,14 @@
 #' This should not typically be called directly
 #'
 #' @param e environment with setup information for popEd
+#' @param eglobal global environment for poped info
 #' @param full setup the full model
 #' @return nothing, called for side effects
-#' @export
 #' @keywords internal
 #' @author Matthew L. Fidler
-.popedSetup <- function(e, full=FALSE) {
-  invisible(.Call(`_babelmixr2_popedSetup`, e, full))
+.popedSetup <- function(e, eglobal, full=FALSE) {
+  popedMultipleEndpointResetTimeIndex()
+  invisible(.Call(`_babelmixr2_popedSetup`, e, eglobal, full))
 }
 #' Solve poped problem for appropriate times with single/multiple endpoint models
 #'
@@ -289,6 +290,7 @@ rxUiGet.popedFgFun <- function(x, ...) {
 attr(rxUiGet.popedFgFun, "desc") <- "PopED parameter model (fg_fun)"
 
 .poped <- new.env(parent=emptyenv())
+.poped$curLib <- NULL
 .poped$s <- NULL
 .poped$modelNumber <- 1L
 .poped$curNumber <- -1L
@@ -418,7 +420,7 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)"
   }
   if (.poped$setup != 1L) {
     rxode2::rxSolveFree()
-    .popedSetup(popedDb$babelmixr2, FALSE)
+    .popedSetup(popedDb$babelmixr2, .poped, FALSE)
     .poped$curNumber <- popedDb$babelmixr2$modelNumber
     .poped$setup <- 1L
     .poped$fullXt <- NULL
@@ -488,7 +490,7 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)"
                                 })))
     .et <- rxode2::etTrans(.dat, .e$modelF)
     .e$dataF <- .et
-    .popedSetup(.e, TRUE)
+    .popedSetup(.e, .poped, TRUE)
     .poped$fullXt <- length(xt)
     .poped$curNumber <- popedDb$babelmixr2$modelNumber
     .poped$setup <- 2L
@@ -555,7 +557,7 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)"
                                 })))
     .et <- rxode2::etTrans(.dat, .e$modelF)
     .e$dataF <- .et
-    .popedSetup(.e, TRUE)
+    .popedSetup(.e, .poped, TRUE)
     .poped$fullXt <- length(xt)
     .poped$curNumber <- popedDb$babelmixr2$modelNumber
     .poped$setup <- 2L
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index abe37bda..1fff9d19 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -94,14 +94,15 @@ BEGIN_RCPP
 END_RCPP
 }
 // popedSetup
-RObject popedSetup(Environment e, bool full);
-RcppExport SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP fullSEXP) {
+RObject popedSetup(Environment e, Environment eglobal, bool full);
+RcppExport SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP eglobalSEXP, SEXP fullSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< Environment >::type e(eSEXP);
+    Rcpp::traits::input_parameter< Environment >::type eglobal(eglobalSEXP);
     Rcpp::traits::input_parameter< bool >::type full(fullSEXP);
-    rcpp_result_gen = Rcpp::wrap(popedSetup(e, full));
+    rcpp_result_gen = Rcpp::wrap(popedSetup(e, eglobal, full));
     return rcpp_result_gen;
 END_RCPP
 }
diff --git a/src/init.c b/src/init.c
index f26569b5..4601685e 100644
--- a/src/init.c
+++ b/src/init.c
@@ -13,7 +13,7 @@ SEXP _babelmixr2_transDv(SEXP, SEXP, SEXP, SEXP, SEXP,
 SEXP _babelmixr2_iniRxodePtrs(SEXP in);
 
 SEXP _babelmixr2_popedFree(void);
-SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP fullSEXP);
+SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP e2SEXP, SEXP fullSEXP);
 SEXP _babelmixr2_popedSolveIdME(SEXP thetaSEXP, SEXP umtSEXP, SEXP mtSEXP, SEXP msSEXP, SEXP nendSEXP, SEXP idSEXP, SEXP totnSEXP);
 SEXP _babelmixr2_popedSolveIdME2(SEXP thetaSEXP, SEXP umtSEXP, SEXP mtSEXP, SEXP msSEXP, SEXP nendSEXP, SEXP idSEXP, SEXP totnSEXP);
 SEXP _babelmixr2_popedGetMultipleEndpointModelingTimes(SEXP, SEXP, SEXP);
@@ -29,7 +29,7 @@ static const R_CallMethodDef CallEntries[] = {
   {"_babelmixr2_popedGetMultipleEndpointModelingTimes",
    (DL_FUNC) &_babelmixr2_popedGetMultipleEndpointModelingTimes, 3},
   {"_babelmixr2_popedFree", (DL_FUNC) &_babelmixr2_popedFree, 0},
-  {"_babelmixr2_popedSetup", (DL_FUNC) &_babelmixr2_popedSetup, 2},
+  {"_babelmixr2_popedSetup", (DL_FUNC) &_babelmixr2_popedSetup, 3},
   {"_babelmixr2_popedSolveIdME", (DL_FUNC) &_babelmixr2_popedSolveIdME, 7},
   {"_babelmixr2_popedSolveIdME2", (DL_FUNC) &_babelmixr2_popedSolveIdME2, 7},
   {"_babelmixr2_iniRxodePtrs", (DL_FUNC) &_babelmixr2_iniRxodePtrs, 1},
diff --git a/src/poped.cpp b/src/poped.cpp
index b3116229..9ca203fe 100644
--- a/src/poped.cpp
+++ b/src/poped.cpp
@@ -23,6 +23,7 @@
 
 timeIndexer globalTimeIndexer;
 
+
 //' @title Get Multiple Endpoint Modeling Times
 //'
 //' @description
@@ -247,6 +248,9 @@ using namespace Rcpp;
 
 #define popedOde(id) ind_solve(rx, id, rxInner.dydt_liblsoda, rxInner.dydt_lsoda_dum, rxInner.jdum_lsoda, rxInner.dydt, rxInner.update_inis, rxInner.global_jt)
 
+Environment _popedE;
+Environment _popedEglobal;
+
 struct rxSolveF {
   //
   // std::string estStr;
@@ -338,12 +342,12 @@ RObject popedFree() {
   return R_NilValue;
 }
 
-Environment _popedE;
 
 //[[Rcpp::export]]
-RObject popedSetup(Environment e, bool full) {
+RObject popedSetup(Environment e, Environment eglobal, bool full) {
   popedFree();
   _popedE=e;
+  _popedEglobal=eglobal;
   List control = e["control"];
   List rxControl = as<List>(e["rxControl"]);
 
@@ -364,7 +368,9 @@ RObject popedSetup(Environment e, bool full) {
   e["paramCache"]=p2;
   e["lid"] = NA_INTEGER;
   List mvp = rxode2::rxModelVars_(model);
-  rxUpdateFuns(as<SEXP>(mvp["trans"]), &rxInner);
+  CharacterVector trans =  mvp["trans"];
+  _popedEglobal["curTrans"] = trans;
+  rxUpdateFuns(as<SEXP>(trans), &rxInner);
 
   // initial value of parameters
   CharacterVector pars = mvp[RxMv_params];
@@ -447,6 +453,7 @@ static inline bool solveCached(NumericVector &theta, int &id) {
 }
 
 void popedSolveFidMat(arma::mat &matMT, NumericVector &theta, int id, int nrow, int nend) {
+  rxUpdateFuns(_popedEglobal["curTrans"], &rxInner);
   // arma::vec ret(retD, nobs, false, true);
   rx_solving_options_ind *ind =  updateParamRetInd(theta, id);
   rx_solving_options *op = getSolvingOptions(rx);
@@ -455,12 +462,12 @@ void popedSolveFidMat(arma::mat &matMT, NumericVector &theta, int id, int nrow,
   int kk, k=0;
   double curT;
   bool isMT = false;
+  double *lhs = getIndLhs(ind);
   for (int j = 0; j < getIndNallTimes(ind); ++j) {
     setIndIdx(ind, j);
     kk = getIndIx(ind, j);
     curT = getTime(kk, ind);
     isMT = getIndEvid(ind, kk) >= 10 && getIndEvid(ind, kk) <= 99;
-    double *lhs = getIndLhs(ind);
     if (isDose(getIndEvid(ind, kk))) {
       rxInner.calc_lhs(id, curT, getOpIndSolve(op, ind, j), lhs);
       continue;
@@ -496,6 +503,7 @@ Rcpp::DataFrame popedSolveIdME(NumericVector &theta,
                                NumericVector &mt, IntegerVector &ms,
                                int nend, int id, int totn) {
   if (solveCached(theta, id)) return(as<Rcpp::DataFrame>(_popedE["s"]));
+  rxUpdateFuns(_popedEglobal["curTrans"], &rxInner);
   NumericVector t(totn);
   arma::vec f(totn);
   arma::vec w(totn);
@@ -508,7 +516,6 @@ Rcpp::DataFrame popedSolveIdME(NumericVector &theta,
     std::fill(curLV.begin(), curLV.end(), 0);
     we[i] = curLV;
   }
-
   popedSolveFidMat(matMT, theta, id, nrow, nend);
   // this gets the information from:
   // - the model time
diff --git a/tests/testthat/test-poped.R b/tests/testthat/test-poped.R
index c4f1596b..07d88e01 100644
--- a/tests/testthat/test-poped.R
+++ b/tests/testthat/test-poped.R
@@ -368,6 +368,76 @@ if (requireNamespace("PopED", quietly=TRUE)) {
 
   })
 
+  test_that("mixing 2 poped databases at the same time", {
+    library(PopED)
+
+    # Define the model
+    f <- function() {
+      ini({
+        tV <- 72.8
+        tKa <- 0.25
+        tCl <- 3.75
+        tF <- fix(0.9)
+        pedCL <- 0.8
+
+        eta.v ~ 0.09
+        eta.ka ~ 0.09
+        eta.cl ~0.25^2
+
+        prop.sd <- fix(sqrt(0.04))
+        add.sd <- fix(sqrt(5e-6))
+
+      })
+      model({
+        V<-tV*exp(eta.v)
+        KA<-tKa*exp(eta.ka) * (pedCL**isPediatric) # add covariate for pediatrics
+        CL<-tCl*exp(eta.cl)
+        Favail <- tF
+
+        N <-  floor(t/TAU)+1
+        y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) *
+          (exp(-CL/V * (t - (N - 1) * TAU)) *
+             (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
+             exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
+
+        y ~ prop(prop.sd) + add(add.sd)
+      })
+    }
+
+    e <- et(c( 1,8,10,240,245))
+
+    babel.db <- nlmixr2(f, e, "poped",
+                        popedControl(m = 2,
+                                     groupsize=20,
+                                     bUseGrouped_xt=TRUE,
+                                     a=list(c(DOSE=20,TAU=24,isPediatric = 0),
+                                            c(DOSE=40, TAU=24,isPediatric = 0))))
+
+
+    # Note, to be able to use the adults FIM to combine with the pediatrics,
+    # both have to have the parameter "pedCL" defined and set notfixed_bpop to 1.
+
+    ## Define pediatric model/design (isPediatric = 1)
+    ## One arm, 4 time points only
+
+    e.ped <- et(c( 1,2,6,240))
+
+    babel.db.ped <- nlmixr2(f, e.ped, "poped",
+                            popedControl(m = 1,
+                                         groupsize=6,
+                                         bUseGrouped_xt=TRUE,
+                                         a=list(c(DOSE=40,TAU=24,isPediatric = 1))))
+
+
+    ##  Create plot of model of adult data without variability
+    vdiffr::expect_doppelganger("poped-adult-first",
+                                plot_model_prediction(babel.db, model_num_points = 300))
+
+    vdiffr::expect_doppelganger("poped-ped-next",
+                                plot_model_prediction(babel.db.ped,
+                                                      model_num_points = 300))
+  })
+
   ##  The tests run interactively runs OK
   ## However, the capture output seems to be interfere with the tests (from PopED)
   # So... these are commented out for now.