diff --git a/inst/include/rxode2parseGetTime.h b/inst/include/rxode2parseGetTime.h
index 0f62221ae..01feba1c2 100644
--- a/inst/include/rxode2parseGetTime.h
+++ b/inst/include/rxode2parseGetTime.h
@@ -38,9 +38,16 @@ static inline double getLag(rx_solving_options_ind *ind, int id, int cmt, double
   }
   double ret = LAG(id, cmt, time);
   if (ISNA(ret)) {
-    op->badSolve=1;
-    if (op->naTime == 0) {
-      op->naTime = 1 + 10*cmt;
+    int newBadSolve = 1;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 1 + 10*cmt;
+#pragma omp atomic write
+      op->naTime = newNaTime;
     }
   }
   return ret;
@@ -51,9 +58,16 @@ static inline double getRate(rx_solving_options_ind *ind, int id, int cmt, doubl
   returnBadTime(t);
   double ret = RATE(id, cmt, dose, t);
   if (ISNA(ret)){
-    op->badSolve=1;
-    if (op->naTime == 0) {
-      op->naTime = 2 + 10*cmt;
+    int newBadSolve = 1;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 2 + 10*cmt;
+#pragma omp atomic write
+      op->naTime = newNaTime;
     }
   }
   return ret;
@@ -65,9 +79,16 @@ static inline double getDur(rx_solving_options_ind *ind, int id, int cmt, double
   if (ISNA(t)) return t;
   double ret = DUR(id, cmt, dose, t);
   if (ISNA(ret)){
-    op->badSolve=1;
-    if (op->naTime == 0) {
-      op->naTime = 3 + 10*cmt;
+    int newBadSolve = 1;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 3 + 10*cmt;
+#pragma omp atomic write
+      op->naTime = newNaTime;
     }
   }
   return ret;
@@ -340,9 +361,16 @@ static inline double handleInfusionItem(int idx, rx_solve *rx, rx_solving_option
     ind->idx = oIdx;
     if (ISNA(f)){
       rx_solving_options *op = &op_global;
-      op->badSolve=1;
-      if (op->naTime == 0) {
-        op->naTime = 4 + 10*ind->cmt;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      int curNaTime;
+#pragma omp atomic read
+      curNaTime = op->naTime;
+      if (curNaTime == 0) {
+        int newNaTime = 4 + 10*ind->cmt;
+#pragma omp atomic write
+        op->naTime = newNaTime;
       }
     }
     double durOld = (getAllTimes(ind, ind->idose[infEidx]) -
diff --git a/inst/include/rxode2parseHandleEvid.h b/inst/include/rxode2parseHandleEvid.h
index f90ed9143..312c29a17 100644
--- a/inst/include/rxode2parseHandleEvid.h
+++ b/inst/include/rxode2parseHandleEvid.h
@@ -285,8 +285,17 @@ static inline double getAmt(rx_solving_options_ind *ind, int id, int cmt,
   double ret = AMT(id, cmt, dose, t, y);
   if (ISNA(ret)){
     rx_solving_options *op = &op_global;
-    op->badSolve=1;
-    op->naTime = 5 + 10*cmt;
+    int newBadSolve = 1;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 5 + 10*cmt;
+#pragma omp atomic write
+      op->naTime = newNaTime;
+    }
   }
   return ret;
 }
@@ -309,20 +318,22 @@ static inline int pushIgnoredDose(int doseIdx, rx_solving_options_ind *ind) {
   }
   if (ind->ignoredDosesN[0]+1 >= ind->ignoredDosesAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->ignoredDoses, (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = 1;
-        // return 0;
+        allocFailed = 1;
       } else {
         ind->ignoredDoses = tmpI;
         ind->ignoredDosesAllocN[0] = (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE);
         re = 1;
       }
     }
-    if (op->badSolve) {
-      return 0; // don't continue if we have a bad solve.
+    if (allocFailed) {
+#pragma omp atomic write
+      op->badSolve = 1;
+      return 0;
     }
   }
   ind->ignoredDoses[ind->ignoredDosesN[0]] = doseIdx;
@@ -334,19 +345,21 @@ static inline int pushPendingDose(int doseIdx, rx_solving_options_ind *ind) {
   int re = 0;
   if (ind->pendingDosesN[0]+1 >= ind->pendingDosesAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->pendingDoses, (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = 1;
-        //return 0;
+        allocFailed = 1;
       } else {
         ind->pendingDoses = tmpI;
         ind->pendingDosesAllocN[0] = (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE);
         re = 1;
       }
     }
-    if (op->badSolve == 1)  {
+    if (allocFailed) {
+#pragma omp atomic write
+      op->badSolve = 1;
       return 0;
     }
   }
@@ -361,32 +374,29 @@ static inline int pushDosingEvent(double time, double amt, int evid,
   int re = 0;
   if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;  // 0=ok, 1=partial alloc, -1=first alloc failed
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = -1;
-        // return 0;
+        allocFailed = -1;
       } else {
         ind->extraDoseTimeIdx = tmpI;
 
         tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
         if (tmpI == NULL) {
-          op->badSolve = 1;
-          // return 1;
+          allocFailed = 1;
         } else {
           ind->extraDoseEvid = tmpI;
           double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
           if (tmpD == NULL) {
-            op->badSolve = 1;
-            //return 1;
+            allocFailed = 1;
           } else {
             ind->extraDoseTime = tmpD;
 
             tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
             if (tmpD == NULL) {
-              op->badSolve = 1;
-              //return 1;
+              allocFailed = 1;
             } else {
               ind->extraDoseDose = tmpD;
 
@@ -397,11 +407,11 @@ static inline int pushDosingEvent(double time, double amt, int evid,
         }
       }
     }
-    if (op->badSolve == 1) {
-      return 1;
-    } else if (op->badSolve == -1) {
-      op->badSolve = 1; // set to bad solve.
-      return 0;
+    if (allocFailed != 0) {
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return (allocFailed == 1) ? 1 : 0;
     }
   }
   ind->extraDoseTimeIdx[ind->extraDoseN[0]] = ind->extraDoseN[0];
@@ -419,33 +429,30 @@ static inline int pushUniqueDosingEvent(double time, double amt, int evid,
   int re = 0;
   if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;  // 0=ok, 1=partial alloc, -1=first alloc failed
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = -1;
-        // return 0;
+        allocFailed = -1;
       }  else {
         ind->extraDoseTimeIdx = tmpI;
 
         tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
         if (tmpI == NULL) {
-          op->badSolve = 1;
-          // return 1;
+          allocFailed = 1;
         } else {
           ind->extraDoseEvid = tmpI;
 
           double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
           if (tmpD == NULL) {
-            op->badSolve = 1;
-            //return 1;
+            allocFailed = 1;
           } else {
             ind->extraDoseTime = tmpD;
 
             tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
             if (tmpD == NULL) {
-              op->badSolve = 1;
-              // return 1;
+              allocFailed = 1;
             } else {
               ind->extraDoseDose = tmpD;
 
@@ -455,11 +462,11 @@ static inline int pushUniqueDosingEvent(double time, double amt, int evid,
         }
       }
     }
-    if (op->badSolve == 1) {
-      return 1; // don't continue if we have a bad solve.
-    } else if (op->badSolve == -1) {
-      op->badSolve = 1; // set to bad solve.
-      return 0;
+    if (allocFailed != 0) {
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return (allocFailed == 1) ? 1 : 0;
     }
     re = 1;
   }
diff --git a/src/dop853.c b/src/dop853.c
index 4716f5d89..f735d933e 100644
--- a/src/dop853.c
+++ b/src/dop853.c
@@ -15,54 +15,9 @@
 
 #define _(String) (String)
 
-static long int     nfcn, nstep, naccpt, nrejct;
-static double       hout, xold, xout;
-static int nrds, *indir;
-static double       *yy1, *k1, *k2, *k3, *k4, *k5, *k6, *k7, *k8, *k9, *k10;
-static double       *rcont1, *rcont2, *rcont3, *rcont4;
-static double       *rcont5, *rcont6, *rcont7, *rcont8;
-
-
-long int nfcnRead (void)
-{
-  return nfcn;
-
-} /* nfcnRead */
-
-
-long int nstepRead (void)
-{
-  return nstep;
-
-} /* stepRead */
-
-
-long int naccptRead (void)
-{
-  return naccpt;
-
-} /* naccptRead */
-
-
-long int nrejctRead (void)
-{
-  return nrejct;
-
-} /* nrejct */
-
-
-double hRead (void)
-{
-  return hout;
-
-} /* hRead */
-
-
-double xRead (void)
-{
-  return xout;
-
-} /* xRead */
+/* All former static variables are now in dop853_ctx_t (see dop853.h).
+   This makes the solver reentrant and thread-safe: each call to dop853()
+   uses its own local context on the stack. */
 
 
 static double sign (double a, double b)
@@ -162,8 +117,36 @@ static double hinit (int *nptr, FcnEqDiff fcn, double x, double* y,
 } /* hinit */
 
 
+/* dense output function */
+static double contd8 (dop853_ctx_t *ctx, int ii, double x)
+{
+  int i;
+  double       s, s1;
+
+  i = INT_MAX;
+
+  if (!ctx->indir)
+    i = ii;
+  else
+    i = ctx->indir[ii];
+
+  if (i == INT_MAX)
+    {
+      Rprintf (_("no dense output available for %uth component"), ii);
+      return 0.0;
+    }
+
+  s = (x - ctx->xold) / ctx->hout;
+  s1 = 1.0 - s;
+
+  return ctx->rcont1[i]+s*(ctx->rcont2[i]+s1*(ctx->rcont3[i]+s*(ctx->rcont4[i]+s1*(ctx->rcont5[i]+
+							       s*(ctx->rcont6[i]+s1*(ctx->rcont7[i]+s*ctx->rcont8[i]))))));
+
+} /* contd8 */
+
+
 /* core integrator */
-static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
+static int dopcor (dop853_ctx_t *ctx, int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                    double hmax, double h, double* rtoler, double* atoler,
                    int itoler, FILE* fileout, SolTrait solout, int iout,
                    long int nmax, double uround, int meth, long int nstiff, double safe,
@@ -191,6 +174,16 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
   double   d61, d66, d67, d68, d69, d610, d611, d612, d613, d614, d615, d616;
   double   d71, d76, d77, d78, d79, d710, d711, d712, d713, d714, d715, d716;
 
+  /* local aliases for context arrays */
+  double *yy1 = ctx->yy1, *k1 = ctx->k1, *k2 = ctx->k2, *k3 = ctx->k3;
+  double *k4 = ctx->k4, *k5 = ctx->k5, *k6 = ctx->k6, *k7 = ctx->k7;
+  double *k8 = ctx->k8, *k9 = ctx->k9, *k10 = ctx->k10;
+  double *rcont1 = ctx->rcont1, *rcont2 = ctx->rcont2;
+  double *rcont3 = ctx->rcont3, *rcont4 = ctx->rcont4;
+  double *rcont5 = ctx->rcont5, *rcont6 = ctx->rcont6;
+  double *rcont7 = ctx->rcont7, *rcont8 = ctx->rcont8;
+  int nrds = ctx->nrds;
+
   /* initialisations */
   switch (meth)
     {
@@ -384,16 +377,16 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
   iord = 8;
   if (h == 0.0)
     h = hinit (nptr, fcn, x, y, posneg, k1, k2, k3, iord, hmax, atoler, rtoler, itoler);
-  nfcn += 2;
+  ctx->nfcn += 2;
   reject = 0;
-  xold = x;
+  ctx->xold = x;
 
   if (iout)
     {
       irtrn = 1;
-      hout = 1.0;
-      xout = x;
-      solout (naccpt+1, xold, x, y, nptr, &irtrn);
+      ctx->hout = 1.0;
+      ctx->xout = x;
+      solout (ctx->naccpt+1, ctx->xold, x, y, nptr, &irtrn);
       if (irtrn < 0)
         {
           /* if (fileout) */
@@ -405,12 +398,12 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
   /* basic integration step */
   while (1)
     {
-      if (nstep > nmax)
+      if (ctx->nstep > nmax)
         {
           /* if (fileout) */
 	  Rprintf (_("exit of dop853 at x = %.16e, more than nmax = %li are needed\n"), x, nmax);
-          xout = x;
-          hout = h;
+          ctx->xout = x;
+          ctx->hout = h;
           return -2;
         }
 
@@ -418,8 +411,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
         {
           /* if (fileout) */
 	  Rprintf (_("exit of dop853 at x = %.16e, step size too small h = %.16e\n"), x, h);
-          xout = x;
-          hout = h;
+          ctx->xout = x;
+          ctx->hout = h;
           return -3;
         }
 
@@ -429,7 +422,7 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           last = 1;
         }
 
-      nstep++;
+      ctx->nstep++;
 
       /* the twelve stages */
       for (i = 0; i < n; i++)
@@ -472,7 +465,7 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                              a127*k7[i] + a128*k8[i] + a129*k9[i] +
                              a1210*k10[i] + a1211*k2[i]);
       fcn (nptr, xph, yy1, k3);
-      nfcn += 11;
+      ctx->nfcn += 11;
       for (i = 0; i < n; i++)
         {
           k4[i] = b1*k1[i] + b6*k6[i] + b7*k7[i] + b8*k8[i] + b9*k9[i] +
@@ -525,12 +518,12 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           /* step accepted */
 
           facold = max_d (err, 1.0E-4);
-          naccpt++;
+          ctx->naccpt++;
           fcn (nptr, xph, k5, k4);
-          nfcn++;
+          ctx->nfcn++;
 
           /* stiffness detection */
-          if (!(naccpt % nstiff) || (iasti > 0))
+          if (!(ctx->naccpt % nstiff) || (iasti > 0))
             {
               stnum = 0.0;
               stden = 0.0;
@@ -550,8 +543,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                   if (iasti == 15)
                     {
 		      Rprintf (_("the problem seems to become stiff at x = %.16e\n"), x);
-		      xout = x;
-		      hout = h;
+		      ctx->xout = x;
+		      ctx->hout = h;
 		      return -4;
                     }
                 }
@@ -621,7 +614,7 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                                      a169*k9[i] + a1613*k4[i] + a1614*k10[i] +
                                      a1615*k2[i]);
               fcn (nptr, x+c16*h, yy1, k3);
-              nfcn += 3;
+              ctx->nfcn += 3;
 
               /* final preparation */
               if (nrds == n)
@@ -653,14 +646,14 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
 
           memcpy (k1, k4, n * sizeof(double));
           memcpy (y, k5, n * sizeof(double));
-          xold = x;
+          ctx->xold = x;
           x = xph;
 
           if (iout)
             {
-              hout = h;
-              xout = x;
-              solout (naccpt+1, xold, x, y, nptr, &irtrn);
+              ctx->hout = h;
+              ctx->xout = x;
+              solout (ctx->naccpt+1, ctx->xold, x, y, nptr, &irtrn);
               if (irtrn < 0)
                 {
                   /* if (fileout) */
@@ -672,8 +665,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           /* normal exit */
           if (last)
             {
-              hout=hnew;
-              xout = x;
+              ctx->hout=hnew;
+              ctx->xout = x;
               return 1;
             }
 
@@ -689,8 +682,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           /* step rejected */
           hnew = h / min_d (facc1, fac11/safe);
           reject = 1;
-          if (naccpt >= 1)
-            nrejct=nrejct + 1;
+          if (ctx->naccpt >= 1)
+            ctx->nrejct=ctx->nrejct + 1;
           last = 0;
         }
 
@@ -707,14 +700,22 @@ int dop853
  double safe, double fac1, double fac2, double beta, double hmax, double h,
  long int nmax, int meth, long int nstiff, int nrdens, int* icont, int licont)
 {
+  dop853_ctx_t ctx;
   int          arret, idid;
   int i;
   int n = nptr[0];
 
-  /* initialisations */
-  nfcn = nstep = naccpt = nrejct = arret = 0;
-  rcont1 = rcont2 = rcont3 = rcont4 = rcont5 = rcont6 = rcont7 = rcont8 = NULL;
-  indir = NULL;
+  /* initialise context */
+  ctx.nfcn = ctx.nstep = ctx.naccpt = ctx.nrejct = 0;
+  ctx.hout = ctx.xold = ctx.xout = 0.0;
+  ctx.nrds = 0;
+  ctx.indir = NULL;
+  ctx.rcont1 = ctx.rcont2 = ctx.rcont3 = ctx.rcont4 = NULL;
+  ctx.rcont5 = ctx.rcont6 = ctx.rcont7 = ctx.rcont8 = NULL;
+  ctx.yy1 = ctx.k1 = ctx.k2 = ctx.k3 = ctx.k4 = ctx.k5 = NULL;
+  ctx.k6 = ctx.k7 = ctx.k8 = ctx.k9 = ctx.k10 = NULL;
+
+  arret = 0;
 
   /* n, the dimension of the system */
   if (n == INT_MAX)
@@ -768,19 +769,19 @@ int dop853
   else if (nrdens)
     {
       /* is there enough memory to allocate rcont12345678&indir ? */
-      rcont1 = R_Calloc(nrdens,double);
-      rcont2 = R_Calloc(nrdens,double);
-      rcont3 = R_Calloc(nrdens,double);
-      rcont4 = R_Calloc(nrdens,double);
-      rcont5 = R_Calloc(nrdens,double);
-      rcont6 = R_Calloc(nrdens,double);
-      rcont7 = R_Calloc(nrdens,double);
-      rcont8 = R_Calloc(nrdens,double);
+      ctx.rcont1 = R_Calloc(nrdens,double);
+      ctx.rcont2 = R_Calloc(nrdens,double);
+      ctx.rcont3 = R_Calloc(nrdens,double);
+      ctx.rcont4 = R_Calloc(nrdens,double);
+      ctx.rcont5 = R_Calloc(nrdens,double);
+      ctx.rcont6 = R_Calloc(nrdens,double);
+      ctx.rcont7 = R_Calloc(nrdens,double);
+      ctx.rcont8 = R_Calloc(nrdens,double);
       if (nrdens < n)
-        indir = R_Calloc(n,int);
+        ctx.indir = R_Calloc(n,int);
 
-      if (!rcont1 || !rcont2 || !rcont3 || !rcont4 || !rcont5 ||
-          !rcont6 || !rcont7 || !rcont8 || (!indir && (nrdens < n)))
+      if (!ctx.rcont1 || !ctx.rcont2 || !ctx.rcont3 || !ctx.rcont4 || !ctx.rcont5 ||
+          !ctx.rcont6 || !ctx.rcont7 || !ctx.rcont8 || (!ctx.indir && (nrdens < n)))
         {
           /* if (fileout) */
 	  Rprintf ( _("not enough free memory for rcont12345678&indir\n"));
@@ -792,7 +793,7 @@ int dop853
         {
           if (icont)
             Rprintf ( _("warning : when nrdens = n there is no need allocating memory for icont\n"));
-          nrds = n;
+          ctx.nrds = n;
         }
       else if (licont < nrdens)
         {
@@ -804,11 +805,11 @@ int dop853
         {
           if ((iout < 2) && fileout)
             fprintf (fileout, "Warning : put iout = 2 for dense output\n");
-          nrds = nrdens;
+          ctx.nrds = nrdens;
           for (i = 0; i < n; i++)
-            indir[i] = INT_MAX;
+            ctx.indir[i] = INT_MAX;
           for (i = 0; i < nrdens; i++)
-            indir[icont[i]] = i;
+            ctx.indir[icont[i]] = i;
         }
     }
 
@@ -855,19 +856,20 @@ int dop853
     hmax = xend - x;
 
   /* is there enough free memory for the method ? */
-  yy1 =R_Calloc(n,double);
-  k1 = R_Calloc(n,double);
-  k2 = R_Calloc(n,double);
-  k3 = R_Calloc(n,double);
-  k4 = R_Calloc(n,double);
-  k5 = R_Calloc(n,double);
-  k6 = R_Calloc(n,double);
-  k7 = R_Calloc(n,double);
-  k8 = R_Calloc(n,double);
-  k9 = R_Calloc(n,double);
-  k10 = R_Calloc(n,double);
-
-  if (!yy1 || !k1 || !k2 || !k3 || !k4 || !k5 || !k6 || !k7 || !k8 || !k9 || !k10)
+  ctx.yy1 =R_Calloc(n,double);
+  ctx.k1 = R_Calloc(n,double);
+  ctx.k2 = R_Calloc(n,double);
+  ctx.k3 = R_Calloc(n,double);
+  ctx.k4 = R_Calloc(n,double);
+  ctx.k5 = R_Calloc(n,double);
+  ctx.k6 = R_Calloc(n,double);
+  ctx.k7 = R_Calloc(n,double);
+  ctx.k8 = R_Calloc(n,double);
+  ctx.k9 = R_Calloc(n,double);
+  ctx.k10 = R_Calloc(n,double);
+
+  if (!ctx.yy1 || !ctx.k1 || !ctx.k2 || !ctx.k3 || !ctx.k4 || !ctx.k5 ||
+      !ctx.k6 || !ctx.k7 || !ctx.k8 || !ctx.k9 || !ctx.k10)
     {
       if (fileout)
         fprintf (fileout, "Not enough free memory for the method\n");
@@ -877,107 +879,79 @@ int dop853
   /* when a failure has occured, we return -1 */
   if (arret)
     {
-      if (k10)
-        R_Free (k10);
-      if (k9)
-        R_Free (k9);
-      if (k8)
-        R_Free (k8);
-      if (k7)
-        R_Free (k7);
-      if (k6)
-        R_Free (k6);
-      if (k5)
-        R_Free (k5);
-      if (k4)
-        R_Free (k4);
-      if (k3)
-        R_Free (k3);
-      if (k2)
-        R_Free (k2);
-      if (k1)
-        R_Free (k1);
-      if (yy1)
-        R_Free (yy1);
-      if (indir)
-        R_Free (indir);
-      if (rcont8)
-        R_Free (rcont8);
-      if (rcont7)
-        R_Free (rcont7);
-      if (rcont6)
-        R_Free (rcont6);
-      if (rcont5)
-        R_Free (rcont5);
-      if (rcont4)
-        R_Free (rcont4);
-      if (rcont3)
-        R_Free (rcont3);
-      if (rcont2)
-        R_Free (rcont2);
-      if (rcont1)
-        R_Free (rcont1);
+      if (ctx.k10)
+        R_Free (ctx.k10);
+      if (ctx.k9)
+        R_Free (ctx.k9);
+      if (ctx.k8)
+        R_Free (ctx.k8);
+      if (ctx.k7)
+        R_Free (ctx.k7);
+      if (ctx.k6)
+        R_Free (ctx.k6);
+      if (ctx.k5)
+        R_Free (ctx.k5);
+      if (ctx.k4)
+        R_Free (ctx.k4);
+      if (ctx.k3)
+        R_Free (ctx.k3);
+      if (ctx.k2)
+        R_Free (ctx.k2);
+      if (ctx.k1)
+        R_Free (ctx.k1);
+      if (ctx.yy1)
+        R_Free (ctx.yy1);
+      if (ctx.indir)
+        R_Free (ctx.indir);
+      if (ctx.rcont8)
+        R_Free (ctx.rcont8);
+      if (ctx.rcont7)
+        R_Free (ctx.rcont7);
+      if (ctx.rcont6)
+        R_Free (ctx.rcont6);
+      if (ctx.rcont5)
+        R_Free (ctx.rcont5);
+      if (ctx.rcont4)
+        R_Free (ctx.rcont4);
+      if (ctx.rcont3)
+        R_Free (ctx.rcont3);
+      if (ctx.rcont2)
+        R_Free (ctx.rcont2);
+      if (ctx.rcont1)
+        R_Free (ctx.rcont1);
 
       return -1;
     }
   else
     {
-      idid = dopcor (nptr, fcn, x, y, xend, hmax, h, rtoler, atoler, itoler, fileout,
+      idid = dopcor (&ctx, nptr, fcn, x, y, xend, hmax, h, rtoler, atoler, itoler, fileout,
                      solout, iout, nmax, uround, meth, nstiff, safe, beta, fac1, fac2, icont);
-      R_Free (k10);
-      R_Free (k9);
-      R_Free (k8);
-      R_Free (k7);
-      R_Free (k6);
-      R_Free (k5);    /* reverse order R_Freeing too increase chances */
-      R_Free (k4);    /* of efficient dynamic memory managing       */
-      R_Free (k3);
-      R_Free (k2);
-      R_Free (k1);
-      R_Free (yy1);
-      if (indir)
-        R_Free (indir);
-      if (rcont8)
+      R_Free (ctx.k10);
+      R_Free (ctx.k9);
+      R_Free (ctx.k8);
+      R_Free (ctx.k7);
+      R_Free (ctx.k6);
+      R_Free (ctx.k5);    /* reverse order R_Freeing too increase chances */
+      R_Free (ctx.k4);    /* of efficient dynamic memory managing       */
+      R_Free (ctx.k3);
+      R_Free (ctx.k2);
+      R_Free (ctx.k1);
+      R_Free (ctx.yy1);
+      if (ctx.indir)
+        R_Free (ctx.indir);
+      if (ctx.rcont8)
         {
-          R_Free (rcont8);
-          R_Free (rcont7);
-          R_Free (rcont6);
-          R_Free (rcont5);
-          R_Free (rcont4);
-          R_Free (rcont3);
-          R_Free (rcont2);
-          R_Free (rcont1);
+          R_Free (ctx.rcont8);
+          R_Free (ctx.rcont7);
+          R_Free (ctx.rcont6);
+          R_Free (ctx.rcont5);
+          R_Free (ctx.rcont4);
+          R_Free (ctx.rcont3);
+          R_Free (ctx.rcont2);
+          R_Free (ctx.rcont1);
         }
 
       return idid;
     }
 
 } /* dop853 */
-
-
-/* dense output function */
-double contd8 (int ii, double x)
-{
-  int i;
-  double       s, s1;
-
-  i = INT_MAX;
-
-  if (!indir)
-    i = ii;
-  else
-    i = indir[ii];
-
-  if (i == INT_MAX)
-    {
-      Rprintf (_("no dense output available for %uth component"), ii);
-      return 0.0;
-    }
-
-  s = (x - xold) / hout;
-  s1 = 1.0 - s;
-
-  return rcont1[i]+s*(rcont2[i]+s1*(rcont3[i]+s*(rcont4[i]+s1*(rcont5[i]+
-							       s*(rcont6[i]+s1*(rcont7[i]+s*rcont8[i]))))));
-
-} /* contd8 */
diff --git a/src/dop853.h b/src/dop853.h
index cbb50223b..ac09f7e66 100644
--- a/src/dop853.h
+++ b/src/dop853.h
@@ -23,8 +23,8 @@ Version of Mai 2, 1994.
 Remarks about the C version : this version allocates memory by itself, the
 iwork array (among the initial FORTRAN parameters) has been splitted into
 independant initial parameters, the statistical variables and last step size
-and x have been encapsulated in the module and are now accessible through
-dedicated functions; the variable names have been kept to maintain a kind
+and x have been encapsulated in a dop853_ctx_t context struct for thread-safe
+reentrant usage; the variable names have been kept to maintain a kind
 of reading compatibility between the C and FORTRAN codes; adaptation made by
 J.Colinge (COLINGE@DIVSUN.UNIGE.CH).
 
@@ -74,10 +74,11 @@ solout   A pointer to the output function called during integration.
 	 Continuous output : during the calls to solout, a continuous solution
 	 for the interval (xold,x) is available through the function
 
-	   contd8(i,s)
+	   contd8(ctx, i, s)
 
 	 which provides an approximation to the i-th component of the solution
-	 at the point s (s must lie in the interval (xold,x)).
+	 at the point s (s must lie in the interval (xold,x)), where ctx is
+	 the dop853_ctx_t context pointer.
 
 iout     Switch for calling solout :
 	   iout=0 : no call,
@@ -145,7 +146,7 @@ Memory requirements
 OUTPUT PARAMETERS
 -----------------
 
-y       numerical solution at x=xRead() (see below).
+y       numerical solution at x=ctx->xout (see dop853_ctx_t below).
 
 dopri5 returns the following values
 
@@ -154,21 +155,22 @@ dopri5 returns the following values
 	-1 : input is not consistent,
 	-2 : larger nmax is needed,
 	-3 : step size becomes too small,
-	-4 : the problem is probably stff (interrupted).
+	-4 : the problem is probably stiff (interrupted).
 
 
-Several functions provide access to different values :
+After a successful call to dop853, the following fields of the dop853_ctx_t
+context struct provide solver statistics and output values :
 
-xRead   x value for which the solution has been computed (x=xend after
+ctx->xout    x value for which the solution has been computed (x=xend after
 	successful return).
 
-hRead   Predicted step size of the last accepted step (useful for a subsequent
-	call to dop853).
+ctx->hout    Predicted step size of the last accepted step (useful for a
+	subsequent call to dop853).
 
-nstepRead   Number of used steps.
-naccptRead  Number of accepted steps.
-nrejctRead  Number of rejected steps.
-nfcnRead    Number of function calls.
+ctx->nstep   Number of used steps.
+ctx->naccpt  Number of accepted steps.
+ctx->nrejct  Number of rejected steps.
+ctx->nfcn    Number of function calls.
 
 
 */
@@ -180,6 +182,17 @@ nfcnRead    Number of function calls.
 typedef void (*FcnEqDiff)(int *nptr, double x, double *y, double *f);
 typedef void (*SolTrait)(long int nr, double xold, double x, double* y, int *nptr, int* irtrn);
 
+/* Context struct holding all solver state, enabling thread-safe reentrant calls.
+   Previously these were file-scope static variables in dop853.c. */
+typedef struct {
+  long int     nfcn, nstep, naccpt, nrejct;
+  double       hout, xold, xout;
+  int          nrds, *indir;
+  double       *yy1, *k1, *k2, *k3, *k4, *k5, *k6, *k7, *k8, *k9, *k10;
+  double       *rcont1, *rcont2, *rcont3, *rcont4;
+  double       *rcont5, *rcont6, *rcont7, *rcont8;
+} dop853_ctx_t;
+
 
 extern int dop853
  (int *nptr,      /* dimension of the system <= INT_MAX-1*/
@@ -207,15 +220,3 @@ extern int dop853
   int* icont, /* indexes of components for which dense output is required, >= nrdens */
   int licont  /* declared length of icon */
  );
-
-extern double contd8
- (int ii,     /* index of desired component */
-  double x             /* approximation at x */
- );
-
-extern long int nfcnRead (void);   /* encapsulation of statistical data */
-extern long int nstepRead (void);
-extern long int naccptRead (void);
-extern long int nrejctRead (void);
-extern double hRead (void);
-extern double xRead (void);
diff --git a/src/handle_evid.cpp b/src/handle_evid.cpp
index 02e760aa3..585f48b90 100644
--- a/src/handle_evid.cpp
+++ b/src/handle_evid.cpp
@@ -31,24 +31,39 @@ extern "C" void handleTlast(double *time, rx_solving_options_ind *ind) {
 }
 
 // Linear compartment models/functions
+// Note: Rf_errorcall is not thread-safe and cannot be called from
+// within OpenMP parallel regions.  Return NA_REAL and set badSolve
+// instead so the error is handled after the parallel region completes.
 extern "C" double _getDur(int l, rx_solving_options_ind *ind, int backward, unsigned int *p) {
   double dose = getDoseNumber(ind, l);
   if (backward==1 && l != 0){
     if (l <= 0) {
-      (Rf_errorcall)(R_NilValue, _("could not find a start to the infusion #1"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     p[0] = l-1;
     while (p[0] > 0 && getDoseNumber(ind, p[0]) != -dose){
       p[0]--;
     }
     if (getDoseNumber(ind, p[0]) != -dose){
-      (Rf_errorcall)(R_NilValue, _("could not find a start to the infusion #2"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     return getAllTimes(ind, ind->idose[l]) - getAllTimes(ind, ind->idose[p[0]]);
   } else {
     if (l >= ind->ndoses) {
       if (backward==2) return(NA_REAL);
-      (Rf_errorcall)(R_NilValue, _("could not find an end to the infusion"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     p[0] = l+1;
     while (p[0] < ind->ndoses && getDoseNumber(ind, p[0]) != -dose){
@@ -56,7 +71,11 @@ extern "C" double _getDur(int l, rx_solving_options_ind *ind, int backward, unsi
     }
     if (getDoseNumber(ind, p[0]) != -dose){
       if (backward==2) return(NA_REAL);
-      (Rf_errorcall)(R_NilValue, _("could not find an end to the infusion"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     return getAllTimes(ind, ind->idose[p[0]]) - getAllTimes(ind, ind->idose[l]);
   }
diff --git a/src/linCmt.cpp b/src/linCmt.cpp
index cd0336d49..67eacfddf 100644
--- a/src/linCmt.cpp
+++ b/src/linCmt.cpp
@@ -40,36 +40,40 @@ extern "C" void ensureLinCmtA(int nCores) {
 }
 
 // Global linear compartment B model object
-// Since this cannot be threaded, this is not a vector
-// object.  This is created once to reduce memory allocation
-// and deallocation time.
-stan::math::linCmtStan __linCmtB(0, 0, 0, true, 0, 0);
-// Maximum size can be 2*ncmt + 1;
-double __linCmtBdata[14];
-int __linCmtBnumSens=0;
+// Refactored to per-thread vector for thread safety, matching linCmtA pattern.
+typedef struct {
+  stan::math::linCmtStan lc;
+  double data[14];
+  int numSens;
+  Eigen::Matrix<double, Eigen::Dynamic, 1> fx;
+  Eigen::Matrix<double, Eigen::Dynamic, 1> yp;
+  Eigen::Matrix<double, Eigen::Dynamic, 2> g;
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> J;
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Js;
+  Eigen::Matrix<double, Eigen::Dynamic, 1> Jg;
+} linB_t;
+
+std::vector<linB_t> __linCmtB;
+
+extern "C" void ensureLinCmtB(int nCores) {
+  if ((int)__linCmtB.size() < nCores) {
+    __linCmtB.resize(nCores);
+  }
+}
 
 #define linCmtBaddrTheta 0
 #define linCmtBaddrThetaSens 1
-static inline double * getLinCmtDoubleAddr(int type) {
+static inline double * getLinCmtDoubleAddr(linB_t &lcb, int type) {
   switch (type) {
   case linCmtBaddrTheta: // max 7
-    return __linCmtBdata;
+    return lcb.data;
   case linCmtBaddrThetaSens:  // max 7
-    return &__linCmtBdata[7];
+    return &lcb.data[7];
   // note fx needs cannot be a Map for use in stan :(
   }
   return NULL;
 }
 
-Eigen::Matrix<double, Eigen::Dynamic, 1> __linCmtBfx;
-
-Eigen::Matrix<double, Eigen::Dynamic, 1> __linCmtByp;
-Eigen::Matrix<double, Eigen::Dynamic, 2> __linCmtBg;
-
-Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> __linCmtBJ;
-Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> __linCmtBJs;
-Eigen::Matrix<double, Eigen::Dynamic, 1> __linCmtBJg;
-
 // [[Rcpp::export]]
 RObject linCmtModelDouble(double dt,
                           double p1, double v1, double p2,
@@ -383,16 +387,16 @@ extern "C" double linCmtA(rx_solve *rx, int id,
 }
 
 extern "C" double linCmtScaleInitPar(int which) {
-  return __linCmtB.initPar(which);
+  return __linCmtB[0].lc.initPar(which);
 }
 
 extern "C" double linCmtScaleInitN() {
-  Eigen::Matrix<double, Eigen::Dynamic, 1> theta = __linCmtB.initPar();
+  Eigen::Matrix<double, Eigen::Dynamic, 1> theta = __linCmtB[0].lc.initPar();
   return theta.size();
 }
 
 extern "C" int linCmtZeroJac(int i) {
-  return __linCmtB.parDepV1(i);
+  return __linCmtB[0].lc.parDepV1(i);
 }
 
 
@@ -481,14 +485,15 @@ extern "C" double linCmtB(rx_solve *rx, int id,
                           double p4, double p5,
                           // Oral parameters
                           double ka) {
-#define fx        __linCmtBfx
-#define Jg        __linCmtBJg
-#define lc        __linCmtB
-#define AlastA    __linCmtBAlastA
-#define J         __linCmtBJ
-#define Js        __linCmtBJs
-#define yp        __linCmtByp
-#define g         __linCmtBg
+  // Per-thread linCmtB state (matching linCmtA pattern for thread safety)
+  linB_t &lcb = __linCmtB[rx_get_thread(__linCmtB.size())];
+#define fx        lcb.fx
+#define Jg        lcb.Jg
+#define lc        lcb.lc
+#define J         lcb.J
+#define Js        lcb.Js
+#define yp        lcb.yp
+#define g         lcb.g
   rx_solving_options_ind *ind = &(rx->subjects[id]);
   rx_solving_options *op = rx->op;
   int idx = ind->idx;
@@ -523,8 +528,8 @@ extern "C" double linCmtB(rx_solve *rx, int id,
     // NA fill and resize
     J = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Constant(ncmt + oral0, npars, NA_REAL);
 
-    __linCmtBnumSens = lc.numSens();
-    Js = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(ncmt+oral0, __linCmtBnumSens);//(ncmt + oral0, 2*ncmt + oral0);
+    lcb.numSens = lc.numSens();
+    Js = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(ncmt+oral0, lcb.numSens);//(ncmt + oral0, 2*ncmt + oral0);
     // thetaSens.resize(numSens);
 
     // AlastA.resize(ncmt + oral0);
@@ -542,7 +547,7 @@ extern "C" double linCmtB(rx_solve *rx, int id,
   lc.setId(id);
 
   Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, 1> >
-    theta(getLinCmtDoubleAddr(linCmtBaddrTheta), lc.getNpars());
+    theta(getLinCmtDoubleAddr(lcb, linCmtBaddrTheta), lc.getNpars());
 
   int sw = ncmt + 10*oral0;
   switch (sw) {
@@ -555,7 +560,7 @@ extern "C" double linCmtB(rx_solve *rx, int id,
   }
 
   Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, 1> >
-    thetaSens(getLinCmtDoubleAddr(linCmtBaddrThetaSens), __linCmtBnumSens);
+    thetaSens(getLinCmtDoubleAddr(lcb, linCmtBaddrThetaSens), lcb.numSens);
 
   lc.sensTheta(theta, thetaSens, rx->sensType == 3, rx->linCmtScale);
   if (ind->linSS == linCmtSsInf) {
@@ -675,13 +680,13 @@ extern "C" double linCmtB(rx_solve *rx, int id,
       }
     }
   }
-  lc.getJacCp(__linCmtBJ, fx, theta, Jg);
+  lc.getJacCp(J, fx, theta, Jg);
   return lc.adjustF(fx, theta, ind->linCmtHV);
 #undef fx
 #undef J
 #undef Jg
 #undef lc
-#undef theta
-#undef AlastA
+#undef Js
 #undef yp
+#undef g
 }
diff --git a/src/par_solve.cpp b/src/par_solve.cpp
index 2c692605c..318eb2958 100644
--- a/src/par_solve.cpp
+++ b/src/par_solve.cpp
@@ -37,6 +37,7 @@ extern "C" {
 #define badSolveExit(i) for (int j = op->neq*(ind->n_all_times); j--;){ \
     ind->solve[j] = NA_REAL;                                           \
   }                                                                     \
+  _Pragma("omp atomic write")                                           \
   op->badSolve = 1;                                                     \
   i = ind->n_all_times-1; // Get out of here!
 // Yay easy parallel support
@@ -2229,6 +2230,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
   double *yp;
   inits = op->inits;
   int idid = 0;
+  int localBadSolve = 0;
   ind = &(rx->subjects[neq[1]]);
   if (!iniSubject(neq[1], 0, ind, op, rx, u_inis)) return;
   nx = ind->n_all_times;
@@ -2249,6 +2251,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
         *rc = -1000;
         // Bad Solve => NA
         badSolveExit(i);
+        localBadSolve = 1;
       } else {
         preSolve(op, ind, xoutp, xout, yp);
         idid = indLin(solveid, op, xoutp, yp, xout, ind->InfusionRate, ind->on,
@@ -2258,7 +2261,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
       }
     }
     ind->_newind = 2;
-    if (!op->badSolve){
+    if (!localBadSolve){
       ind->idx = i;
       if (getEvid(ind, ind->ix[i]) == 3) {
         handleEvid3(ind, op, rx, neq, &xp, &xout,  yp, &idid, u_inis);
@@ -2301,7 +2304,8 @@ extern "C" void par_indLin(rx_solve *rx){
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   // FIXME parallel
   uint32_t seed0 = getRxSeed1(1);
   for (int solveid = 0; solveid < nsim*nsub; solveid++){
@@ -2356,6 +2360,7 @@ extern "C" void ind_liblsoda0(rx_solve *rx, rx_solving_options *op, struct lsoda
   double *yp;
   int neqOde = *neq - op->numLin - op->numLinSens;
   inits = op->inits;
+  int localBadSolve = 0;
   struct lsoda_context_t * ctx = lsoda_create_ctx();
   ctx->function = (_lsoda_f)dydt_liblsoda;
   ctx->data = neq;
@@ -2386,45 +2391,51 @@ extern "C" void ind_liblsoda0(rx_solve *rx, rx_solving_options *op, struct lsoda
         *rc = -1000;
         // Bad Solve => NA
         badSolveExit(i);
+        localBadSolve = 1;
       } else {
         // REprintf("xp: %f xout: %f\n", xp, xout);
         if (handleExtraDose(neq, BadDose, InfusionRate, ind->dose, yp, xout,
                             xp, ind->id, &i, nx, &(ctx->state), op, ind, u_inis, ctx)) {
-          if (!isSameTime(ind->extraDoseNewXout, xp)) {
+          if (!localBadSolve && !isSameTime(ind->extraDoseNewXout, xp)) {
             preSolve(op, ind, xp, ind->extraDoseNewXout, yp);
             lsoda(ctx, yp, &xp, ind->extraDoseNewXout);
             copyLinCmt(neq, ind, op, yp);
             postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+            if (*rc < 0) localBadSolve = 1;
           }
-          int idx = ind->idx;
-          int ixds= ind->ixds;
-          int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
-          ind->idx = -1-trueIdx;
-          handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
-                      BadDose, InfusionRate, ind->dose, yp, xout, neq[1], ind);
-          ctx->state = 1;
-          ind->idx = idx;
-          ind->ixds = ixds;
-          ind->idxExtra++;
-          if (!isSameTime(xout, ind->extraDoseNewXout)) {
-            preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
-            lsoda(ctx,yp, &ind->extraDoseNewXout, xout);
-            copyLinCmt(neq, ind, op, yp);
-            postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+          if (!localBadSolve) {
+            int idx = ind->idx;
+            int ixds= ind->ixds;
+            int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
+            ind->idx = -1-trueIdx;
+            handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
+                        BadDose, InfusionRate, ind->dose, yp, xout, neq[1], ind);
+            ctx->state = 1;
+            ind->idx = idx;
+            ind->ixds = ixds;
+            ind->idxExtra++;
+            if (!isSameTime(xout, ind->extraDoseNewXout)) {
+              preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
+              lsoda(ctx,yp, &ind->extraDoseNewXout, xout);
+              copyLinCmt(neq, ind, op, yp);
+              postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+              if (*rc < 0) localBadSolve = 1;
+            }
+            xp =  ind->extraDoseNewXout;
           }
-          xp =  ind->extraDoseNewXout;
         }
-        if (!isSameTime(xout, xp)) {
+        if (!localBadSolve && !isSameTime(xout, xp)) {
           preSolve(op, ind, xp, xout, yp);
           lsoda(ctx, yp, &xp, xout);
           copyLinCmt(neq, ind, op, yp);
           postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+          if (*rc < 0) localBadSolve = 1;
         }
         xp = xout;
       }
     }
     ind->_newind = 2;
-    if (!op->badSolve){
+    if (!localBadSolve){
       ind->idx = i;
       if (getEvid(ind, ind->ix[i]) == 3) {
         handleEvid3(ind, op, rx, neq, &xp, &xout,  yp, &(ctx->state), u_inis);
@@ -2494,7 +2505,8 @@ extern "C" void par_linCmt(rx_solve *rx) {
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   uint32_t seed0 = getRxSeed1(cores);
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(cores)
@@ -2571,7 +2583,8 @@ extern "C" void par_liblsodaR(rx_solve *rx) {
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   uint32_t seed0 = getRxSeed1(cores);
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(cores)
@@ -2646,7 +2659,8 @@ extern "C" void par_liblsoda(rx_solve *rx){
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   uint32_t seed0 = getRxSeed1(cores);
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(op->cores)
@@ -2816,6 +2830,7 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
 
   double xp = getAllTimes(ind, 0);
   double xout;
+  int localBadSolve = 0;
 
   if (!iniSubject(neq[1], 0, ind, op, rx, u_inis)) return;
   ind->solvedIdx = 0;
@@ -2829,10 +2844,11 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
         ind->rc[0] = -1000;
         // Bad Solve => NA
         badSolveExit(i);
+        localBadSolve = 1;
       } else {
         if (handleExtraDose(neq, ind->BadDose, ind->InfusionRate, ind->dose, yp, xout,
                             xp, ind->id, &i, ind->n_all_times, &istate, op, ind, u_inis, ctx)) {
-          if (!isSameTime(ind->extraDoseNewXout, xp)) {
+          if (!localBadSolve && !isSameTime(ind->extraDoseNewXout, xp)) {
             preSolve(op, ind, xp, ind->extraDoseNewXout, yp);
             neq[0] = op->neq - op->numLin - op->numLinSens;
             F77_CALL(dlsoda)(dydt_lsoda, neq, yp, &xp, &ind->extraDoseNewXout, &gitol, &(op->RTOL), &(op->ATOL), &gitask,
@@ -2840,29 +2856,33 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
             neq[0] = op->neq;
             copyLinCmt(neq, ind, op, yp);
             postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+            if (*(ind->rc) < 0) localBadSolve = 1;
           }
-          int idx = ind->idx;
-          int ixds = ind->ixds;
-          int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
-          ind->idx = -1-trueIdx;
-          handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
-                      ind->BadDose, ind->InfusionRate, ind->dose, yp, xout, neq[1], ind);
-          istate = 1;
-          ind->ixds = ixds; // This is a fake dose, real dose stays in place
-          ind->idx = idx;
-          ind->idxExtra++;
-          if (!isSameTime(xout, ind->extraDoseNewXout)) {
-            preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
-            neq[0] = op->neq - op->numLin - op->numLinSens;
-            F77_CALL(dlsoda)(dydt_lsoda, neq, yp, &ind->extraDoseNewXout, &xout, &gitol, &(op->RTOL), &(op->ATOL), &gitask,
-                             &istate, &giopt, rwork, &lrw, iwork, &liw, jdum, &jt);
-            neq[0] = op->neq;
-            copyLinCmt(neq, ind, op, yp);
-            postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+          if (!localBadSolve) {
+            int idx = ind->idx;
+            int ixds = ind->ixds;
+            int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
+            ind->idx = -1-trueIdx;
+            handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
+                        ind->BadDose, ind->InfusionRate, ind->dose, yp, xout, neq[1], ind);
+            istate = 1;
+            ind->ixds = ixds; // This is a fake dose, real dose stays in place
+            ind->idx = idx;
+            ind->idxExtra++;
+            if (!isSameTime(xout, ind->extraDoseNewXout)) {
+              preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
+              neq[0] = op->neq - op->numLin - op->numLinSens;
+              F77_CALL(dlsoda)(dydt_lsoda, neq, yp, &ind->extraDoseNewXout, &xout, &gitol, &(op->RTOL), &(op->ATOL), &gitask,
+                               &istate, &giopt, rwork, &lrw, iwork, &liw, jdum, &jt);
+              neq[0] = op->neq;
+              copyLinCmt(neq, ind, op, yp);
+              postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+              if (*(ind->rc) < 0) localBadSolve = 1;
+            }
+            xp =  ind->extraDoseNewXout;
           }
-          xp =  ind->extraDoseNewXout;
         }
-        if (!isSameTime(xout, xp)) {
+        if (!localBadSolve && !isSameTime(xout, xp)) {
           preSolve(op, ind, xp, xout, yp);
           neq[0] = op->neq - op->numLin - op->numLinSens;
           F77_CALL(dlsoda)(dydt_lsoda, neq, yp,
@@ -2874,13 +2894,14 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
           neq[0] = op->neq;
           copyLinCmt(neq, ind, op, yp);
           postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+          if (*(ind->rc) < 0) localBadSolve = 1;
         }
         xp = xout;
         //dadt_counter = 0;
       }
     }
     ind->_newind = 2;
-    if (!op->badSolve){
+    if (!localBadSolve){
       ind->idx = i;
       if (getEvid(ind, ind->ix[i]) == 3) {
         handleEvid3(ind, op, rx, neq, &xp, &xout,  yp, &(istate), u_inis);
diff --git a/src/par_solve.h b/src/par_solve.h
index dad810589..a9363a2fb 100644
--- a/src/par_solve.h
+++ b/src/par_solve.h
@@ -78,7 +78,10 @@ extern "C" {
 			ind->solved = -1;
 		}
     sortInd(ind);
-    if (op->badSolve) return 0;
+    // Note: previously checked op->badSolve here, but that is a shared
+    // global flag — in parallel mode another thread's failed solve would
+    // prevent THIS individual from initializing.  The bad-solve state is
+    // now handled per-individual via localBadSolve / *rc in the caller.
 		ind->ixds=ind->idx=0;
     if (ncmt) ind->pendingDosesN[0] = 0;
 		return 1;
diff --git a/src/rxData.cpp b/src/rxData.cpp
index 2b9aabc88..1c3f7d84d 100644
--- a/src/rxData.cpp
+++ b/src/rxData.cpp
@@ -43,6 +43,7 @@ using namespace arma;
 
 extern "C" void seedEng(int ncores);
 extern "C" void ensureLinCmtA(int nCores);
+extern "C" void ensureLinCmtB(int nCores);
 
 #include "cbindThetaOmega.h"
 #include "../inst/include/rxode2parseHandleEvid.h"
@@ -4994,6 +4995,7 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
     }
     seedEng((int)(op->cores));
     ensureLinCmtA((int)op->cores);
+    ensureLinCmtB((int)op->cores);
 
     // Now set up events and parameters
     RObject par0 = params;
diff --git a/src/rxomp.h b/src/rxomp.h
index 308f3630a..6d67bb410 100644
--- a/src/rxomp.h
+++ b/src/rxomp.h
@@ -8,8 +8,10 @@
 static inline int rx_get_thread(int mx) {
   int tn = omp_get_thread_num();
   if (tn < 0) return 0;
-  if (tn <= mx) return tn;
-  return 0;
+  if (tn < mx) return tn;
+  // Clamp to last valid index rather than collapsing to 0,
+  // which would cause two threads to share the same slot
+  return (mx > 0) ? mx - 1 : 0;
 }
 
 #else