diff --git a/inst/include/rxode2parseGetTime.h b/inst/include/rxode2parseGetTime.h
index 2c3af735..aa4a0923 100644
--- a/inst/include/rxode2parseGetTime.h
+++ b/inst/include/rxode2parseGetTime.h
@@ -7,11 +7,11 @@ extern "C" {
 
 #include "rxode2parse.h"
 
-extern t_F AMT;
-extern t_LAG LAG;
-extern t_RATE RATE;
-extern t_DUR DUR;
-extern t_calc_mtime calc_mtime;
+  extern t_F AMT;
+  extern t_LAG LAG;
+  extern t_RATE RATE;
+  extern t_DUR DUR;
+  extern t_calc_mtime calc_mtime;
 
 #ifndef __DOINIT__
 
@@ -30,432 +30,375 @@ extern t_calc_mtime calc_mtime;
   }
 
 
-static inline double getLag(rx_solving_options_ind *ind, int id, int cmt, double time) {
-  rx_solving_options *op = &op_global;
-  returnBadTime(time);
-  if (ind->wh0 == EVID0_SS0 || ind->wh0 == EVID0_SS20) {
-    return time;
-  }
-  double ret = LAG(id, cmt, time);
-  if (ISNA(ret)) {
-    op->badSolve=1;
-    op->naTime = 1;
+  static inline double getLag(rx_solving_options_ind *ind, int id, int cmt, double time) {
+    rx_solving_options *op = &op_global;
+    returnBadTime(time);
+    if (ind->wh0 == EVID0_SS0 || ind->wh0 == EVID0_SS20) {
+      return time;
+    }
+    double ret = LAG(id, cmt, time);
+    if (ISNA(ret)) {
+      op->badSolve=1;
+      op->naTime = 1;
+    }
+    return ret;
   }
-  return ret;
-}
 
-static inline double getRate(rx_solving_options_ind *ind, int id, int cmt, double dose, double t){
-  rx_solving_options *op = &op_global;
-  returnBadTime(t);
-  double ret = RATE(id, cmt, dose, t);
-  if (ISNA(ret)){
-    op->badSolve=1;
-    op->naTime = 1;
+  static inline double getRate(rx_solving_options_ind *ind, int id, int cmt, double dose, double t){
+    rx_solving_options *op = &op_global;
+    returnBadTime(t);
+    double ret = RATE(id, cmt, dose, t);
+    if (ISNA(ret)){
+      op->badSolve=1;
+      op->naTime = 1;
+    }
+    return ret;
   }
-  return ret;
-}
 
-static inline double getDur(rx_solving_options_ind *ind, int id, int cmt, double dose, double t){
-  rx_solving_options *op = &op_global;
-  returnBadTime(t);
-  if (ISNA(t)) return t;
-  double ret = DUR(id, cmt, dose, t);
-  if (ISNA(ret)){
-    op->badSolve=1;
-    op->naTime = 1;
+  static inline double getDur(rx_solving_options_ind *ind, int id, int cmt, double dose, double t){
+    rx_solving_options *op = &op_global;
+    returnBadTime(t);
+    if (ISNA(t)) return t;
+    double ret = DUR(id, cmt, dose, t);
+    if (ISNA(ret)){
+      op->badSolve=1;
+      op->naTime = 1;
+    }
+    return ret;
   }
-  return ret;
-}
 
 
-static inline int isEvidType(int evid, int type) {
-  int wh, cmt, wh100, whI, wh0;
-  getWh(evid, &wh, &cmt, &wh100, &whI, &wh0);
-  return (whI == type);
-}
+  static inline int isEvidType(int evid, int type) {
+    int wh, cmt, wh100, whI, wh0;
+    getWh(evid, &wh, &cmt, &wh100, &whI, &wh0);
+    return (whI == type);
+  }
 
 #define isEvidModeledDurationStart(evid) isEvidType(evid, EVIDF_MODEL_DUR_ON)
 #define isEvidModeledDurationStop(evid) isEvidType(evid, EVIDF_MODEL_DUR_OFF)
 #define isEvidModeledRateStart(evid) isEvidType(evid, EVIDF_MODEL_RATE_ON)
 #define isEvidModeledRateStop(evid) isEvidType(evid, EVIDF_MODEL_RATE_OFF)
 
-static inline void updateDur(int idx, rx_solving_options_ind *ind, double *yp){
-  double t = getAllTimes(ind, idx);
-  double dur, amt;
-  // The duration and f cannot depend on state values
-  int oldIdx = ind->idx;
-  ind->idx = idx;
-  amt  = getAmt(ind, ind->id, ind->cmt, getDose(ind, idx), t, yp);
-  dur  = getDur(ind, ind->id, ind->cmt, amt, t);
-  ind->idx = oldIdx;
-  if (dur > 0) {
-    setDoseP1(ind, idx, -amt/dur);
-    setAllTimesP1(ind, idx, t+dur);
-  } else {
-    rx_solve *rx = &rx_global;
-    rx_solving_options *op = &op_global;
-    if (ind->cmt < op->neq){
-      if (rx->needSort & 4){
-        if (!(ind->err & 16)){
-          ind->err += 16;
-        }
-        return;
-      } else {
-        if (!(ind->err & 32)){
-          ind->err += 32;
+  static inline void updateDur(int idx, rx_solving_options_ind *ind, double *yp){
+    double t = getAllTimes(ind, idx);
+    double dur, amt;
+    // The duration and f cannot depend on state values
+    int oldIdx = ind->idx;
+    ind->idx = idx;
+    amt  = getAmt(ind, ind->id, ind->cmt, getDose(ind, idx), t, yp);
+    dur  = getDur(ind, ind->id, ind->cmt, amt, t);
+    ind->idx = oldIdx;
+    if (dur > 0) {
+      setDoseP1(ind, idx, -amt/dur);
+      setAllTimesP1(ind, idx, t+dur);
+    } else {
+      rx_solve *rx = &rx_global;
+      rx_solving_options *op = &op_global;
+      if (ind->cmt < op->neq){
+        if (rx->needSort & 4){
+          if (!(ind->err & 16)){
+            ind->err += 16;
+          }
+          return;
+        } else {
+          if (!(ind->err & 32)){
+            ind->err += 32;
+          }
+          return;
         }
-        return;
       }
     }
   }
-}
 
-static inline void updateRate(int idx, rx_solving_options_ind *ind, double *yp) {
-  double t = getAllTimes(ind, idx);
-  int oldIdx = ind->idx;
-  ind->idx=idx;
-  double dur, rate, amt;
-  amt  = getAmt(ind, ind->id, ind->cmt, getDose(ind,idx), t, yp);
-  rate  = getRate(ind, ind->id, ind->cmt, amt, t);
-  if (rate > 0){
-    dur = amt/rate; // mg/hr
-    setDoseP1(ind, idx, -rate);
-    setAllTimesP1(ind, idx, t+dur);
-    ind->idx=oldIdx;
-  } else {
-    rx_solve *rx;
-    rx = &rx_global;
-    rx_solving_options *op = &op_global;
-    if (ind->cmt < op->neq){
-      if (rx->needSort & 8){
-        if (!(ind->err & 2)){
-          ind->err += 2;
-        }
-        return;
-      } else {
-        // FIXME don't error out with linear compartmental model
-        if (!(ind->err & 4)){
-          ind->err += 4;
+  static inline void updateRate(int idx, rx_solving_options_ind *ind, double *yp) {
+    double t = getAllTimes(ind, idx);
+    int oldIdx = ind->idx;
+    ind->idx=idx;
+    double dur, rate, amt;
+    amt  = getAmt(ind, ind->id, ind->cmt, getDose(ind,idx), t, yp);
+    rate  = getRate(ind, ind->id, ind->cmt, amt, t);
+    if (rate > 0){
+      dur = amt/rate; // mg/hr
+      setDoseP1(ind, idx, -rate);
+      setAllTimesP1(ind, idx, t+dur);
+      ind->idx=oldIdx;
+    } else {
+      rx_solve *rx;
+      rx = &rx_global;
+      rx_solving_options *op = &op_global;
+      if (ind->cmt < op->neq){
+        if (rx->needSort & 8){
+          if (!(ind->err & 2)){
+            ind->err += 2;
+          }
+          return;
+        } else {
+          // FIXME don't error out with linear compartmental model
+          if (!(ind->err & 4)){
+            ind->err += 4;
+          }
+          return;
         }
-        return;
       }
     }
+    ind->idx=oldIdx;
   }
-  ind->idx=oldIdx;
-}
 
-static inline void handleTurnOffModeledDuration(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
-  if (idx > 0){
-    if (!isEvidModeledDurationStart(getEvidM1(ind, idx))) {
-      if (!(ind->err & 64)){
-        ind->err += 64;
+  static inline void handleTurnOffModeledDuration(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
+    if (idx > 0){
+      if (!isEvidModeledDurationStart(getEvidM1(ind, idx))) {
+        if (!(ind->err & 64)){
+          ind->err += 64;
+        }
+        return;
       }
-      return;
-    }
-  } else {
-    if (!(ind->err & 128)){
-      ind->err += 128;
-    }
-    return;
-  }
-}
-
-static inline void handleTurnOnModeledDuration(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
-  // This calculates the rate and the duration and then assigns it to the next record
-  if (idx >= ind->n_all_times){
-    // error: Last record, can't be used.
-    if (!(ind->err & 256)){
-      ind->err += 256;
-    }
-    return;
-  } else {
-    if (!isEvidModeledDurationStop(getEvidP1(ind, idx))) {
-      if (!(ind->err & 512)){
-        ind->err += 512;
+    } else {
+      if (!(ind->err & 128)){
+        ind->err += 128;
       }
       return;
     }
-    updateDur(idx, ind, rx->ypNA);
   }
-}
 
-static inline void handleTurnOffModeledRate(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
-  if (idx > 0){
-    if (!isEvidModeledRateStart(getEvidM1(ind, idx))) {
-      if (!(ind->err & 1024)){
-        ind->err += 1024;
+  static inline void handleTurnOnModeledDuration(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
+    // This calculates the rate and the duration and then assigns it to the next record
+    if (idx >= ind->n_all_times){
+      // error: Last record, can't be used.
+      if (!(ind->err & 256)){
+        ind->err += 256;
       }
       return;
-    }
-  } else {
-    if (!(ind->err & 2048)){
-      ind->err += 2048;
-    }
-    return;
-  }
-}
-
-static inline void handleTurnOnModeledRate(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
-  // This calculates the rate and the duration and then assigns it to the next record
-  if (idx >= ind->n_all_times){
-    // error: Last record, can't be used.
-    if (!(ind->err & 4096)){
-      ind->err += 4096;
-    }
-    /* Rf_errorcall(R_NilValue, "Data Error 9\n"); */
-    return;
-  } else {
-    if (!isEvidModeledRateStop(getEvidP1(ind, idx))) {
-      if (!(ind->err & 8192)){
-        ind->err += 8192;
+    } else {
+      if (!isEvidModeledDurationStop(getEvidP1(ind, idx))) {
+        if (!(ind->err & 512)){
+          ind->err += 512;
+        }
+        return;
       }
-      return;
+      updateDur(idx, ind, rx->ypNA);
     }
-    setAllTimesP1(ind, idx, getAllTimes(ind, idx));
-    updateRate(idx, ind, rx->ypNA);
   }
-}
 
-static inline int handleInfusionStartRm(int *startIdx, int *endIdx,
-                                        double *amt, int *idx,
-                                        rx_solve *rx, rx_solving_options *op,
-                                        rx_solving_options_ind *ind) {
-  if (ind->wh0 == EVID0_INFRM) {
-    // This is a possible removal event.  Look at the next duration
-    *startIdx = *endIdx+1;
-    for (*endIdx = *startIdx; *endIdx < ind->ndoses; ++(*endIdx)) {
-      if (getEvid(ind, ind->idose[*startIdx]) == getEvid(ind, ind->idose[*endIdx])) break;
-      if (*endIdx == ind->ndoses-1) {
-        //REprintf("curEvid@infrm: %d\n", curEvid);
-        if (!(ind->err & 32768)){
-          ind->err += 32768;
+  static inline void handleTurnOffModeledRate(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
+    if (idx > 0){
+      if (!isEvidModeledRateStart(getEvidM1(ind, idx))) {
+        if (!(ind->err & 1024)){
+          ind->err += 1024;
         }
-        return 1;
+        return;
+      }
+    } else {
+      if (!(ind->err & 2048)){
+        ind->err += 2048;
       }
+      return;
     }
-    return 1;
   }
-  return 0;
-}
 
-static inline int handleInfusionStartDefault(int *startIdx, int *endIdx,
-                                             double *amt, int *idx,
-                                             rx_solve *rx, rx_solving_options *op,
-                                             rx_solving_options_ind *ind) {
-  // This finds the duration based on the end of infusion
-  int curEvid = getEvid(ind, ind->idose[*endIdx]);
-  int jj = 0;
-  for (*startIdx = 0; *startIdx < ind->ndoses; (*startIdx)++) {
-    if (getEvid(ind, ind->idose[*startIdx]) == curEvid &&
-        getDose(ind, ind->idose[*startIdx]) == -(*amt)) {
-      // This will look after the last known infusion stopping point
-      if (jj == 0) {
-        jj = *startIdx;
-      } else {
-        jj++;
+  static inline void handleTurnOnModeledRate(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
+    // This calculates the rate and the duration and then assigns it to the next record
+    if (idx >= ind->n_all_times){
+      // error: Last record, can't be used.
+      if (!(ind->err & 4096)){
+        ind->err += 4096;
       }
-      for (; jj < ind->ndoses; jj++) {
-        if (getEvid(ind, ind->idose[jj]) == curEvid &&
-            getDose(ind, ind->idose[jj]) == *amt) {
-          break;
+      /* Rf_errorcall(R_NilValue, "Data Error 9\n"); */
+      return;
+    } else {
+      if (!isEvidModeledRateStop(getEvidP1(ind, idx))) {
+        if (!(ind->err & 8192)){
+          ind->err += 8192;
         }
+        return;
       }
-      if (jj == *endIdx) break;
+      setAllTimesP1(ind, idx, getAllTimes(ind, idx));
+      updateRate(idx, ind, rx->ypNA);
     }
   }
-  if (*startIdx == ind->ndoses) {
-    //REprintf("cant match: %d\n", curEvid);
-    if (!(ind->err & 32768)){
-      ind->err += 32768;
-    }
-    return 1;
-  }
-  return 1;
-}
 
-static inline void handleInfusionGetStartOfInfusionIndex(int *startIdx, int *endIdx,
-                                                         double *amt, int *idx,
-                                                         rx_solve *rx, rx_solving_options *op,
-                                                         rx_solving_options_ind *ind) {
-  if (*amt > 0) {
-    *startIdx = -1;
-    *endIdx = -1;
-    return;
-  }
-  *endIdx = getDoseNumberFromIndex(ind, *idx);
-  if (*endIdx == -1){
-    if (!(ind->err & 16384)){
-      ind->err += 16384;
+  static inline void handleInfusionGetStartOfInfusionIndex(int *startIdx, int *endIdx,
+                                                           double *amt, int *idx,
+                                                           rx_solve *rx, rx_solving_options *op,
+                                                           rx_solving_options_ind *ind) {
+    if (*amt > 0) {
+      *startIdx = -1;
+      *endIdx = -1;
+      return;
     }
-    return;
-    /* Rf_errorcall(R_NilValue, "Corrupted event table during sort (1)."); */
-  }
-  handleInfusionStartRm(startIdx, endIdx, amt, idx, rx, op, ind) ||
-    handleInfusionStartDefault(startIdx, endIdx, amt, idx, rx, op, ind);
-}
-
-static inline double handleInfusionItem(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
-  if (ind->wh0 == EVID0_RATEADJ) {
-    return getAllTimes(ind, idx);
-  }
-  double amt = getDose(ind, idx);
-  if (amt > 0) {
-    return getLag(ind, ind->id, ind->cmt, getAllTimes(ind, idx));
-  } else if (amt < 0) {
-    int infEidx = getDoseNumberFromIndex(ind, idx);
-    if (infEidx == -1){
+    *endIdx = getDoseNumberFromIndex(ind, *idx);
+    if (*endIdx == -1){
       if (!(ind->err & 16384)){
         ind->err += 16384;
       }
-      return 0.0;
+      return;
       /* Rf_errorcall(R_NilValue, "Corrupted event table during sort (1)."); */
     }
-    int infBidx;
-    handleInfusionGetStartOfInfusionIndex(&infBidx, &infEidx, &amt, &idx, rx, op, ind);
-    if (infBidx == -1) return 0.0;
-    rx_solve *rx = &rx_global;
-    double f = getAmt(ind, ind->id, ind->cmt, 1.0, getAllTimes(ind, ind->idose[infBidx]), rx->ypNA);
-    if (ISNA(f)){
-      rx_solving_options *op = &op_global;
-      op->badSolve=1;
-      op->naTime = 1;
+    handleInfusionStartRm(startIdx, endIdx, amt, idx, rx, op, ind) ||
+      handleInfusionStartDefault(startIdx, endIdx, amt, idx, rx, op, ind);
+  }
+
+  static inline double handleInfusionItem(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
+    if (ind->wh0 == EVID0_RATEADJ) {
+      return getAllTimes(ind, idx);
     }
-    double durOld = (getAllTimes(ind, ind->idose[infEidx]) -
-                     getAllTimes(ind, ind->idose[infBidx]));
-    double dur = f*durOld;
-    // To the correct lag time, we need to make sure the steady state flags are not set...
-    int wh0 = ind->wh0;
-    int wh, cmt, wh100, whI;
-    getWh(getEvid(ind, ind->idose[infBidx]), &wh, &cmt, &wh100, &whI, &(ind->wh0));
-    double tB = getLag(ind, ind->id, ind->cmt, getAllTimes(ind, ind->idose[infBidx]));
-    ind->wh0 = wh0;
-    return tB + dur;
-  } else {
-    /* Rf_errorcall(R_NilValue, "Corrupted events."); */
-    if (!(ind->err & 131072)){
-      ind->err += 131072;
+    double amt = getDose(ind, idx);
+    if (amt > 0) {
+      return getLag(ind, ind->id, ind->cmt, getAllTimes(ind, idx));
+    } else if (amt < 0) {
+      int infEidx = getDoseNumberFromIndex(ind, idx);
+      if (infEidx == -1){
+        if (!(ind->err & 16384)){
+          ind->err += 16384;
+        }
+        return 0.0;
+        /* Rf_errorcall(R_NilValue, "Corrupted event table during sort (1)."); */
+      }
+      int infBidx;
+      handleInfusionGetStartOfInfusionIndex(&infBidx, &infEidx, &amt, &idx, rx, op, ind);
+      if (infBidx == -1) return 0.0;
+      rx_solve *rx = &rx_global;
+      double f = getAmt(ind, ind->id, ind->cmt, 1.0, getAllTimes(ind, ind->idose[infBidx]), rx->ypNA);
+      if (ISNA(f)){
+        rx_solving_options *op = &op_global;
+        op->badSolve=1;
+        op->naTime = 1;
+      }
+      double durOld = (getAllTimes(ind, ind->idose[infEidx]) -
+                       getAllTimes(ind, ind->idose[infBidx]));
+      double dur = f*durOld;
+      // To the correct lag time, we need to make sure the steady state flags are not set...
+      int wh0 = ind->wh0;
+      int wh, cmt, wh100, whI;
+      getWh(getEvid(ind, ind->idose[infBidx]), &wh, &cmt, &wh100, &whI, &(ind->wh0));
+      double tB = getLag(ind, ind->id, ind->cmt, getAllTimes(ind, ind->idose[infBidx]));
+      ind->wh0 = wh0;
+      return tB + dur;
+    } else {
+      /* Rf_errorcall(R_NilValue, "Corrupted events."); */
+      if (!(ind->err & 131072)){
+        ind->err += 131072;
+      }
+      return 0.0;
     }
-    return 0.0;
   }
-}
 
-static inline double getTimeCalculateInfusionTimes(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
-  switch(ind->whI){
-  case EVIDF_MODEL_DUR_OFF:
-    handleTurnOffModeledDuration(idx, rx, op, ind);
-    break;
-  case EVIDF_MODEL_DUR_ON:
-    handleTurnOnModeledDuration(idx, rx, op, ind);
-    break;
-  case EVIDF_MODEL_RATE_OFF:
-    handleTurnOffModeledRate(idx, rx, op, ind);
-    break;
-  case EVIDF_MODEL_RATE_ON:
-    handleTurnOnModeledRate(idx, rx, op, ind);
-    break;
-  case EVIDF_INF_RATE:
-    return handleInfusionItem(idx, rx, op, ind);
-    break;
+  static inline double getTimeCalculateInfusionTimes(int idx, rx_solve *rx, rx_solving_options *op, rx_solving_options_ind *ind) {
+    switch(ind->whI){
+    case EVIDF_MODEL_DUR_OFF:
+      handleTurnOffModeledDuration(idx, rx, op, ind);
+      break;
+    case EVIDF_MODEL_DUR_ON:
+      handleTurnOnModeledDuration(idx, rx, op, ind);
+      break;
+    case EVIDF_MODEL_RATE_OFF:
+      handleTurnOffModeledRate(idx, rx, op, ind);
+      break;
+    case EVIDF_MODEL_RATE_ON:
+      handleTurnOnModeledRate(idx, rx, op, ind);
+      break;
+    case EVIDF_INF_RATE:
+      return handleInfusionItem(idx, rx, op, ind);
+      break;
+    }
+    return getLag(ind, ind->id, ind->cmt, getAllTimes(ind,idx));
   }
-  return getLag(ind, ind->id, ind->cmt, getAllTimes(ind,idx));
-}
 
-static inline double getTime__(int idx, rx_solving_options_ind *ind, int update) {
-  rx_solving_options *op = &op_global;
-  rx_solve *rx = &rx_global;
-  int evid = getEvid(ind, idx);
-  if (evid == 9) return 0.0;
-  if (evid >= 10 && evid <= 99) return ind->mtime[evid-10];
-  if (isObs(evid)) return getAllTimes(ind, idx);
-  getWh(evid, &(ind->wh), &(ind->cmt), &(ind->wh100), &(ind->whI), &(ind->wh0));
-  if (ind->wh0 == EVID0_SSINF ||
-      ind->wh0 == EVID0_SS0 ||
-      ind->wh0 == EVID0_SS20){
-  } else {
-    // yp should be the current solve values
-    //
-    // Before solving the solve will be zero
-    // After solving the yp will contain the solved values
-    //
-    if (update == 0) {
-      if (ind->whI == EVIDF_INF_RATE) {
-        return handleInfusionItem(idx, rx, op, ind);
-      }
+  static inline double getTime__(int idx, rx_solving_options_ind *ind, int update) {
+    rx_solving_options *op = &op_global;
+    rx_solve *rx = &rx_global;
+    int evid = getEvid(ind, idx);
+    if (evid == 9) return 0.0;
+    if (evid >= 10 && evid <= 99) return ind->mtime[evid-10];
+    if (isObs(evid)) return getAllTimes(ind, idx);
+    getWh(evid, &(ind->wh), &(ind->cmt), &(ind->wh100), &(ind->whI), &(ind->wh0));
+    if (ind->wh0 == EVID0_SSINF ||
+        ind->wh0 == EVID0_SS0 ||
+        ind->wh0 == EVID0_SS20){
     } else {
-      return getTimeCalculateInfusionTimes(idx, rx, op, ind);
+      // yp should be the current solve values
+      //
+      // Before solving the solve will be zero
+      // After solving the yp will contain the solved values
+      //
+      if (update == 0) {
+        if (ind->whI == EVIDF_INF_RATE) {
+          return handleInfusionItem(idx, rx, op, ind);
+        }
+      } else {
+        return getTimeCalculateInfusionTimes(idx, rx, op, ind);
+      }
     }
+    return getLag(ind, ind->id, ind->cmt, getAllTimes(ind, idx));
   }
-  return getLag(ind, ind->id, ind->cmt, getAllTimes(ind, idx));
-}
-
-static inline double getTime_(int idx, rx_solving_options_ind *ind) {
-  return getTime__(idx, ind, 0);
-}
 
-static inline int cancelPendingDoses(rx_solving_options_ind *ind, int id) {
-  int re = 0;
-  for (int i = 0; i < ind->pendingDosesN[0]; ++i) {
-    re = pushIgnoredDose(ind->pendingDoses[i], ind) || re;
+  static inline double getTime_(int idx, rx_solving_options_ind *ind) {
+    return getTime__(idx, ind, 0);
   }
-  ind->pendingDosesN[0] = 0;
-  return re;
-}
 
-static inline int cancelInfusionsThatHaveStarted(rx_solving_options_ind *ind, int id, double time) {
-  int re = 0;
-  int wh, cmt, wh100, whI, wh0, evid;
-  double amt, curTime;
-  for (int i = 0; i < ind->extraDoseN[0]; ++i) {
-    int cur = -1-i;
-    evid = getEvid(ind, cur);
-    // REprintf("Extra dose information (%d):\n", cur);
-    // REprintf("\ttime: %f (curTime %f)\n", getAllTimes(ind, cur), time);
-    // REprintf("\tevid: %d ", getEvid(ind, cur));
-    // REprintf("\tdose: %f \n", getDose(ind, cur));
-    // REprintf("================================================================================\n");
-    getWh(evid, &wh, &cmt, &wh100, &whI, &wh0);
-    if (whI != 1 && whI != 2) {
-      // not an infusion
-      continue;
+  static inline int cancelPendingDoses(rx_solving_options_ind *ind, int id) {
+    int re = 0;
+    for (int i = 0; i < ind->pendingDosesN[0]; ++i) {
+      re = pushIgnoredDose(ind->pendingDoses[i], ind) || re;
     }
-    // infusions
-    amt = getDose(ind, cur);
-    curTime = getAllTimes(ind, cur);
-    if (amt > 0) {
-      // an infusion start dose
-      bool ignore = curTime < time;
-      // infusion starts after time;
-      // These should still be pending
-      if (ignore) {
-        // REprintf("ignore!\n");
-        re = pushIgnoredDose(cur, ind) || re;
-      }
-      i++;
-      if (i >= ind->extraDoseN[0]) break;
-      cur = -1-i;
-      double nextAmt = getDose(ind, cur);
+    ind->pendingDosesN[0] = 0;
+    return re;
+  }
+
+  static inline int cancelInfusionsThatHaveStarted(rx_solving_options_ind *ind, int id, double time) {
+    int re = 0;
+    int wh, cmt, wh100, whI, wh0, evid;
+    double amt, curTime;
+    for (int i = 0; i < ind->extraDoseN[0]; ++i) {
+      int cur = -1-i;
       evid = getEvid(ind, cur);
-      int nextWhI =0;
-      getWh(evid, &wh, &cmt, &wh100, &nextWhI, &wh0);
-      // REprintf("pair (%d):\n", cur);
+      // REprintf("Extra dose information (%d):\n", cur);
       // REprintf("\ttime: %f (curTime %f)\n", getAllTimes(ind, cur), time);
-      // REprintf("\tevid: %d ", evid);
-      // REprintf("\tdose: %f \n", nextAmt);
+      // REprintf("\tevid: %d ", getEvid(ind, cur));
+      // REprintf("\tdose: %f \n", getDose(ind, cur));
       // REprintf("================================================================================\n");
-      if (amt == -nextAmt && whI == nextWhI) {
+      getWh(evid, &wh, &cmt, &wh100, &whI, &wh0);
+      if (whI != 1 && whI != 2) {
+        // not an infusion
+        continue;
+      }
+      // infusions
+      amt = getDose(ind, cur);
+      curTime = getAllTimes(ind, cur);
+      if (amt > 0) {
+        // an infusion start dose
+        bool ignore = curTime < time;
+        // infusion starts after time;
+        // These should still be pending
         if (ignore) {
           // REprintf("ignore!\n");
           re = pushIgnoredDose(cur, ind) || re;
         }
-      } else {
-        i--;
+        i++;
+        if (i >= ind->extraDoseN[0]) break;
+        cur = -1-i;
+        double nextAmt = getDose(ind, cur);
+        evid = getEvid(ind, cur);
+        int nextWhI =0;
+        getWh(evid, &wh, &cmt, &wh100, &nextWhI, &wh0);
+        // REprintf("pair (%d):\n", cur);
+        // REprintf("\ttime: %f (curTime %f)\n", getAllTimes(ind, cur), time);
+        // REprintf("\tevid: %d ", evid);
+        // REprintf("\tdose: %f \n", nextAmt);
+        // REprintf("================================================================================\n");
+        if (amt == -nextAmt && whI == nextWhI) {
+          if (ignore) {
+            // REprintf("ignore!\n");
+            re = pushIgnoredDose(cur, ind) || re;
+          }
+        } else {
+          i--;
+        }
+        continue;
       }
-      continue;
+      // turn infusion off. Not sure what to do here. For now keep.
     }
-    // turn infusion off. Not sure what to do here. For now keep.
+    return re;
   }
-  return re;
-}
 
 #undef cancelOrPush
 #undef returnBadTime
diff --git a/inst/include/rxode2parseHandleEvid.h b/inst/include/rxode2parseHandleEvid.h
index 0d5d7064..049ff51a 100644
--- a/inst/include/rxode2parseHandleEvid.h
+++ b/inst/include/rxode2parseHandleEvid.h
@@ -6,7 +6,7 @@ extern "C" {
 #endif
 
 #include "rxode2parse.h"
-//#include "rxThreadData.h"
+  //#include "rxThreadData.h"
 #define rxErrCorruptETSort    1
 #define rxErrRate0            2
 #define rxErrModelRateAbsent  4
@@ -34,40 +34,40 @@ extern "C" {
 
 #if defined(__cplusplus)
 #define FLOOR(x) std::floor(x)
-extern "C" {
+  extern "C" {
 #else
 #define FLOOR(x) floor(x)
 #endif
 
 #ifndef _isrxode2parse_
 
-  int handle_evidL(int evid, double *yp, double xout, rx_solving_options_ind *ind);
-  void handleTlast(double *time, rx_solving_options_ind *ind);
-  double _getDur(int l, rx_solving_options_ind *ind, int backward, unsigned int *p);
+    int handle_evidL(int evid, double *yp, double xout, rx_solving_options_ind *ind);
+    void handleTlast(double *time, rx_solving_options_ind *ind);
+    double _getDur(int l, rx_solving_options_ind *ind, int backward, unsigned int *p);
 
 #else
-  #define _getDur _rxode2parse_getDur
-  extern rx_solving_options op_global;
-  extern rx_solve rx_global;
+#define _getDur _rxode2parse_getDur
+    extern rx_solving_options op_global;
+    extern rx_solve rx_global;
 #endif
 
 #if defined(__cplusplus)
-}
+  }
 #endif
 
 #define EVID_EXTRA_SIZE 10
-// EVID = 0; Observations
-// EVID = 1; is illegal, but converted from NONMEM
-// EVID = 2; Non-observation, possibly covariate
-// EVID = 3; Reset ODE states to zero; Non-observation event
-// EVID = 4; Reset and then dose event;  Illegal
-// EVID = 9; Non-observation event to ini system at time zero; This is to set the INIs at the correct place.
-// EVID = 10-99; mtime events (from ODE system)
-// When EVID > 100
-// EVID: ## # ## ##
-//       c2 I c1 xx
-// c2 = Compartment numbers over 100
-//  I = Infusion Flag/ Special event flag
+  // EVID = 0; Observations
+  // EVID = 1; is illegal, but converted from NONMEM
+  // EVID = 2; Non-observation, possibly covariate
+  // EVID = 3; Reset ODE states to zero; Non-observation event
+  // EVID = 4; Reset and then dose event;  Illegal
+  // EVID = 9; Non-observation event to ini system at time zero; This is to set the INIs at the correct place.
+  // EVID = 10-99; mtime events (from ODE system)
+  // When EVID > 100
+  // EVID: ## # ## ##
+  //       c2 I c1 xx
+  // c2 = Compartment numbers over 100
+  //  I = Infusion Flag/ Special event flag
 #define EVIDF_NORMAL 0
 
 #define EVIDF_INF_RATE 1
@@ -81,28 +81,28 @@ extern "C" {
 
 #define EVIDF_MODEL_RATE_ON  9
 #define EVIDF_MODEL_RATE_OFF 7
-//      0 = no Infusion
-//      1 = Infusion, AMT=rate (mg/hr for instance)
-//      2 = Infusion, duration is fixed
-//      4 = Replacement event
-//      5 = Multiplication event
-//      6 = Turn off modeled duration
-//      7 = Turn off modeled rate compartment
-//      8 = Duration is modeled, AMT=dose; Rate = AMT/(Modeled Duration) NONMEM RATE=-2
-//      9 = Rate is modeled, AMT=dose; Duration = AMT/(Modeled Rate) NONMEM RATE=-1
-// c1 = Compartment numbers below 99
-// xx =  1, regular event (no lag time)
-// xx =  2, An infusion/rate event that doesn't look for start/end of infusion AND does not apply lags
-// xx =  8, possibly turn off steady state infusion with lag time (needed in case spans dur)
-// xx =  9, steady state event SS=1 with lag time
-// xx = 10, steady state event SS=1 (no lag)
-// xx = 19, steady state event at dose time (SS=2) with lag
-// xx = 20, steady state event + last observed info (not lagged)
-// xx = 21, steady state event at dose time (with absorption lag) + last observed info
-// xx = 30, Turn off compartment
-// xx = 40, Steady state constant infusion
-// xx = 50, Phantom event, used for transit compartments
-// Steady state events need a II data item > 0
+  //      0 = no Infusion
+  //      1 = Infusion, AMT=rate (mg/hr for instance)
+  //      2 = Infusion, duration is fixed
+  //      4 = Replacement event
+  //      5 = Multiplication event
+  //      6 = Turn off modeled duration
+  //      7 = Turn off modeled rate compartment
+  //      8 = Duration is modeled, AMT=dose; Rate = AMT/(Modeled Duration) NONMEM RATE=-2
+  //      9 = Rate is modeled, AMT=dose; Duration = AMT/(Modeled Rate) NONMEM RATE=-1
+  // c1 = Compartment numbers below 99
+  // xx =  1, regular event (no lag time)
+  // xx =  2, An infusion/rate event that doesn't look for start/end of infusion AND does not apply lags
+  // xx =  8, possibly turn off steady state infusion with lag time (needed in case spans dur)
+  // xx =  9, steady state event SS=1 with lag time
+  // xx = 10, steady state event SS=1 (no lag)
+  // xx = 19, steady state event at dose time (SS=2) with lag
+  // xx = 20, steady state event + last observed info (not lagged)
+  // xx = 21, steady state event at dose time (with absorption lag) + last observed info
+  // xx = 30, Turn off compartment
+  // xx = 40, Steady state constant infusion
+  // xx = 50, Phantom event, used for transit compartments
+  // Steady state events need a II data item > 0
 #define EVID0_REGULAR  1
 #define EVID0_RATEADJ 2
 #define EVID0_INFRM 8
@@ -115,645 +115,703 @@ extern "C" {
 #define EVID0_PHANTOM 50
 
 
-static inline void getWh(int evid, int *wh, int *cmt, int *wh100, int *whI, int *wh0){
-  *wh = evid;
-  *cmt = 0;
-  *wh100 = FLOOR(*wh/1e5L);
-  *whI   = FLOOR(*wh/1e4L-*wh100*10);
-  *wh    = *wh - *wh100*1e5 - (*whI-1)*1e4;
-  *wh0 = FLOOR((*wh%10000)/100);
-  *cmt = *wh0 - 1 + *wh100*100;
-  *wh0 = evid - *wh100*1e5 - *whI*1e4 - *wh0*100;
-  if (rx_global.linNcmt != 0) {
-    if (rx_global.linKa) {
-      switch (*cmt) {
-      case 0:
-        *cmt = op_global.neq;
-        break;
-      case 1:
-        *cmt = op_global.neq+1;
-        break;
-      case 2:
-        *cmt -= 2;
-        break;
-      }
-    } else {
-      if (*cmt == 0) {
-        *cmt = op_global.neq;
+  static inline void getWh(int evid, int *wh, int *cmt, int *wh100, int *whI, int *wh0){
+    *wh = evid;
+    *cmt = 0;
+    *wh100 = FLOOR(*wh/1e5L);
+    *whI   = FLOOR(*wh/1e4L-*wh100*10);
+    *wh    = *wh - *wh100*1e5 - (*whI-1)*1e4;
+    *wh0 = FLOOR((*wh%10000)/100);
+    *cmt = *wh0 - 1 + *wh100*100;
+    *wh0 = evid - *wh100*1e5 - *whI*1e4 - *wh0*100;
+    if (rx_global.linNcmt != 0) {
+      if (rx_global.linKa) {
+        switch (*cmt) {
+        case 0:
+          *cmt = op_global.neq;
+          break;
+        case 1:
+          *cmt = op_global.neq+1;
+          break;
+        case 2:
+          *cmt -= 2;
+          break;
+        }
       } else {
-        *cmt -= 1;
+        if (*cmt == 0) {
+          *cmt = op_global.neq;
+        } else {
+          *cmt -= 1;
+        }
       }
     }
   }
-}
 
-static inline double getDoseNumber(rx_solving_options_ind *ind, int i) {
-  return getDose(ind, ind->idose[i]);
-}
+  static inline double getDoseNumber(rx_solving_options_ind *ind, int i) {
+    return getDose(ind, ind->idose[i]);
+  }
 
-static inline double getDoseIndex(rx_solving_options_ind *ind, int i) {
-  return (i < 0 ? getDose(ind, i) : getDose(ind, ind->ix[i]));
-}
+  static inline double getDoseIndex(rx_solving_options_ind *ind, int i) {
+    return (i < 0 ? getDose(ind, i) : getDose(ind, ind->ix[i]));
+  }
 
-static inline double getDoseIndexPlus1(rx_solving_options_ind *ind, int i) {
-  return getDoseP1(ind, ind->ix[i]);
-}
+  static inline double getDoseIndexPlus1(rx_solving_options_ind *ind, int i) {
+    return getDoseP1(ind, ind->ix[i]);
+  }
 
-static inline double getIiNumber(rx_solving_options_ind *ind, int i) {
-  //return ind->ii[i];
-  return getIi(ind,ind->idose[i]);
-}
+  static inline double getIiNumber(rx_solving_options_ind *ind, int i) {
+    //return ind->ii[i];
+    return getIi(ind,ind->idose[i]);
+  }
 
-static inline void setDoseNumber(rx_solving_options_ind *ind, int i, int j, double value) {
-  setDoseP1(ind, ind->idose[i] + j, value)
-}
+  static inline void setDoseNumber(rx_solving_options_ind *ind, int i, int j, double value) {
+    setDoseP1(ind, ind->idose[i] + j, value)
+      }
 
-static inline void handleInfusionGetEndOfInfusionIndex(int idx, int *infEixds,
-																											 rx_solve *rx, rx_solving_options *op,
-																											 rx_solving_options_ind *ind) {
-	int curEvid = getEvid(ind, ind->idose[idx]);
-	double curAmt = getDoseNumber(ind, idx);
-	int lastKnownOff = 0;
-	*infEixds = -1;
-	for (int j = 0; j < ind->ndoses; j++) {
-		if (curEvid == getEvid(ind, ind->idose[j]) &&
-				curAmt == getDoseNumber(ind, j)) {
-			// get the first dose combination
-			if (lastKnownOff == 0) {
-				lastKnownOff=j+1;
-			} else {
-				lastKnownOff++;
-			}
-			for (int k = lastKnownOff; k < ind->ndoses; k++) {
-				if (curEvid == getEvid(ind, ind->idose[k]) &&
-						curAmt == -getDoseNumber(ind, k)) {
-					lastKnownOff = k;
-					if (j == idx) {
-						*infEixds = k;
-						// dur = getTime_(ind->idose[infEixds], ind);// -
-						// dur -= getTime_(ind->idose[ind->ixds+2], ind);
-						// dur2 = getIiNumber(ind, ind->ixds) - dur;
-					}
-					k = ind->ndoses;
-				}
-			}
-		}
-		if (*infEixds != -1) break;
-	}
-}
+  static inline void handleInfusionGetEndOfInfusionIndex(int idx, int *infEixds,
+                                                         rx_solve *rx, rx_solving_options *op,
+                                                         rx_solving_options_ind *ind) {
+    int curEvid = getEvid(ind, ind->idose[idx]);
+    double curAmt = getDoseNumber(ind, idx);
+    int lastKnownOff = 0;
+    *infEixds = -1;
+    for (int j = 0; j < ind->ndoses; j++) {
+      if (curEvid == getEvid(ind, ind->idose[j]) &&
+          curAmt == getDoseNumber(ind, j)) {
+        // get the first dose combination
+        if (lastKnownOff == 0) {
+          lastKnownOff=j+1;
+        } else {
+          lastKnownOff++;
+        }
+        for (int k = lastKnownOff; k < ind->ndoses; k++) {
+          if (curEvid == getEvid(ind, ind->idose[k]) &&
+              curAmt == -getDoseNumber(ind, k)) {
+            lastKnownOff = k;
+            if (j == idx) {
+              *infEixds = k;
+              // dur = getTime_(ind->idose[infEixds], ind);// -
+              // dur -= getTime_(ind->idose[ind->ixds+2], ind);
+              // dur2 = getIiNumber(ind, ind->ixds) - dur;
+            }
+            k = ind->ndoses;
+          }
+        }
+      }
+      if (*infEixds != -1) break;
+    }
+  }
 
-static inline int handleTlastInlineUpateDosingInformation(rx_solving_options_ind *ind, double *curDose, double *tinf) {
-  unsigned int p;
-  switch (ind->whI) {
-  case EVIDF_MODEL_RATE_ON: // modeled rate.
-  case EVIDF_MODEL_DUR_ON: // modeled duration.
-    // Rate already calculated and saved in the next dose record
-    // InfusionRate[cmt] -= getDoseIndexPlus1(ind, ind->idx);
-    *tinf = getAllTimesP1(ind, ind->idx) - getAllTimes(ind, ind->idx);
-    return 1;
-    break;
-  case EVIDF_MODEL_RATE_OFF: // End modeled rate
-  case EVIDF_MODEL_DUR_OFF: // end modeled duration
+  static inline int handleInfusionStartRm(int *startIdx, int *endIdx,
+                                          double *amt, int *idx,
+                                          rx_solve *rx, rx_solving_options *op,
+                                          rx_solving_options_ind *ind) {
+    if (ind->wh0 == EVID0_INFRM) {
+      // This is a possible removal event.  Look at the next duration
+      *startIdx = *endIdx+1;
+      for (*endIdx = *startIdx; *endIdx < ind->ndoses; ++(*endIdx)) {
+        if (getEvid(ind, ind->idose[*startIdx]) == getEvid(ind, ind->idose[*endIdx])) break;
+        if (*endIdx == ind->ndoses-1) {
+          //REprintf("curEvid@infrm: %d\n", curEvid);
+          if (!(ind->err & 32768)){
+            ind->err += 32768;
+          }
+          return 1;
+        }
+      }
+      return 1;
+    }
     return 0;
-    break;
-  case EVIDF_INF_DUR:
-  case EVIDF_INF_RATE:
-    if (curDose[0] <= 0) {
+  }
+
+  static inline int handleInfusionStartDefault(int *startIdx, int *endIdx,
+                                               double *amt, int *idx,
+                                               rx_solve *rx, rx_solving_options *op,
+                                               rx_solving_options_ind *ind) {
+    // This finds the duration based on the end of infusion
+    int curEvid = getEvid(ind, ind->idose[*endIdx]);
+    int jj = 0;
+    for (*startIdx = 0; *startIdx < ind->ndoses; (*startIdx)++) {
+      if (getEvid(ind, ind->idose[*startIdx]) == curEvid &&
+          getDose(ind, ind->idose[*startIdx]) == -(*amt)) {
+        // This will look after the last known infusion stopping point
+        if (jj == 0) {
+          jj = *startIdx;
+        } else {
+          jj++;
+        }
+        for (; jj < ind->ndoses; jj++) {
+          if (getEvid(ind, ind->idose[jj]) == curEvid &&
+              getDose(ind, ind->idose[jj]) == *amt) {
+            break;
+          }
+        }
+        if (jj == *endIdx) break;
+      }
+    }
+    if (*startIdx == ind->ndoses) {
+      //REprintf("cant match: %d\n", curEvid);
+      if (!(ind->err & 32768)){
+        ind->err += 32768;
+      }
+      return 1;
+    }
+    return 1;
+  }
+
+  static inline int handleTlastInlineUpateDosingInformation(rx_solving_options_ind *ind, double *curDose, double *tinf) {
+    unsigned int p;
+    switch (ind->whI) {
+    case EVIDF_MODEL_RATE_ON: // modeled rate.
+    case EVIDF_MODEL_DUR_ON: // modeled duration.
+      // Rate already calculated and saved in the next dose record
+      // InfusionRate[cmt] -= getDoseIndexPlus1(ind, ind->idx);
+      *tinf = getAllTimesP1(ind, ind->idx) - getAllTimes(ind, ind->idx);
+      return 1;
+      break;
+    case EVIDF_MODEL_RATE_OFF: // End modeled rate
+    case EVIDF_MODEL_DUR_OFF: // end modeled duration
       return 0;
-    } else {
-      // The amt in rxode2 is the infusion rate, but we need the amt
-      tinf[0] = _getDur(ind->ixds, ind, 2, &p);
-      if (!ISNA(tinf[0])) {
-        curDose[0] = tinf[0] * curDose[0];
-        return 1;
-      } else {
+      break;
+    case EVIDF_INF_DUR:
+    case EVIDF_INF_RATE:
+      if (curDose[0] <= 0) {
         return 0;
+      } else {
+        // The amt in rxode2 is the infusion rate, but we need the amt
+
+        tinf[0] = _getDur(ind->ixds, ind, 2, &p);
+        if (!ISNA(tinf[0])) {
+          curDose[0] = tinf[0] * curDose[0];
+          return 1;
+        } else {
+          return 0;
+        }
       }
+      break;
     }
-    break;
+    return 1;
   }
-  return 1;
-}
 
-static inline void handleTlastInline(double *time, rx_solving_options_ind *ind) {
-  rx_solving_options *op = &op_global;
-  double _time = *time + ind->curShift;
-  int evid = 0;
-  if (ind->idx < 0) {
-    evid = getEvid(ind, ind->idx);
-  } else {
-    evid = getEvid(ind, ind->ix[ind->idx]);
-  }
-  if (op->neq + op->extraCmt != 0 && ind->tlast != _time &&
-      isDose(evid) &&
-      ind->cmt < op->neq + op->extraCmt) {
-    double curDose = getDoseIndex(ind, ind->idx), tinf = NA_REAL;
-    if (handleTlastInlineUpateDosingInformation(ind, &curDose, &tinf) == 0) return;
-    ind->dosenum++;
-    ind->tlast = _time;
-    ind->curDose = curDose;
-    ind->curDoseS[ind->cmt] = ind->curDose;
-    if (ISNA(ind->tfirst)) ind->tfirst = _time;
-    ind->tlastS[ind->cmt] = _time;
-    if (ISNA(ind->tfirstS[ind->cmt])) ind->tfirstS[ind->cmt] = _time;
+  static inline void handleTlastInline(double *time, rx_solving_options_ind *ind) {
+    rx_solving_options *op = &op_global;
+    double _time = *time + ind->curShift;
+    int evid = 0;
+    if (ind->idx < 0) {
+      evid = getEvid(ind, ind->idx);
+    } else {
+      evid = getEvid(ind, ind->ix[ind->idx]);
+    }
+    if (op->neq + op->extraCmt != 0 && ind->tlast != _time &&
+        isDose(evid) &&
+        ind->cmt < op->neq + op->extraCmt) {
+      double curDose = getDoseIndex(ind, ind->idx), tinf = NA_REAL;
+      if (handleTlastInlineUpateDosingInformation(ind, &curDose, &tinf) == 0) return;
+      ind->dosenum++;
+      ind->tlast = _time;
+      ind->curDose = curDose;
+      ind->curDoseS[ind->cmt] = ind->curDose;
+      if (ISNA(ind->tfirst)) ind->tfirst = _time;
+      ind->tlastS[ind->cmt] = _time;
+      if (ISNA(ind->tfirstS[ind->cmt])) ind->tfirstS[ind->cmt] = _time;
+    }
   }
-}
 
-static inline int getDoseNumberFromIndex(rx_solving_options_ind *ind, int idx) {
-  // bisection https://en.wikipedia.org/wiki/Binary_search_algorithm
-  int l = 0, r = ind->ndoses-1, m=0, idose = 0;
-  while(l <= r){
-    m = FLOOR((l+r)/2);
-    idose= ind->idose[m];
-    if (idose < idx) l = m+1;
-    else if (idose > idx) r = m-1;
-    else return m;
+  static inline int getDoseNumberFromIndex(rx_solving_options_ind *ind, int idx) {
+    // bisection https://en.wikipedia.org/wiki/Binary_search_algorithm
+    int l = 0, r = ind->ndoses-1, m=0, idose = 0;
+    while(l <= r){
+      m = FLOOR((l+r)/2);
+      idose= ind->idose[m];
+      if (idose < idx) l = m+1;
+      else if (idose > idx) r = m-1;
+      else return m;
+    }
+    return -1;
   }
-  return -1;
-}
 
 
-static inline int syncIdx(rx_solving_options_ind *ind) {
-  if (ind->idx < 0) return 1; // additional dose; technically the idx doesn't relate to idose/ix
-  if (ind->ix[ind->idx] != ind->idose[ind->ixds]) {
-    // bisection https://en.wikipedia.org/wiki/Binary_search_algorithm
-    int m = getDoseNumberFromIndex(ind, ind->ix[ind->idx]);
-    if (m != -1) {
-      // REprintf("sync to %d ind->ixds from %d to %d #1 (ind->idx: %d)\n", getEvid(ind, ind->idose[m]),
-      //          ind->ixds, m, ind->idx);
-      ind->ixds=m;
-    } else {
-      //262144
-      if (!(ind->err & rxErrSync)){
-        ind->err += rxErrSync;
-      }
-      return 0;
-    }
-    // Need to adjust ixdsr
-    for(int j = ind->ixds; j--;){
-      if (ind->ix[ind->idx] == ind->idose[j]){
-        ind->ixds = j;
-        break;
+  static inline int syncIdx(rx_solving_options_ind *ind) {
+    if (ind->idx < 0) return 1; // additional dose; technically the idx doesn't relate to idose/ix
+    if (ind->ix[ind->idx] != ind->idose[ind->ixds]) {
+      // bisection https://en.wikipedia.org/wiki/Binary_search_algorithm
+      int m = getDoseNumberFromIndex(ind, ind->ix[ind->idx]);
+      if (m != -1) {
+        // REprintf("sync to %d ind->ixds from %d to %d #1 (ind->idx: %d)\n", getEvid(ind, ind->idose[m]),
+        //          ind->ixds, m, ind->idx);
+        ind->ixds=m;
+      } else {
+        //262144
+        if (!(ind->err & rxErrSync)){
+          ind->err += rxErrSync;
+        }
+        return 0;
       }
-    }
-    if (ind->ix[ind->idx] != ind->idose[ind->ixds]){
-      for(int j = ind->ixds+1; j< ind->ndoses; j++){
+      // Need to adjust ixdsr
+      for(int j = ind->ixds; j--;){
         if (ind->ix[ind->idx] == ind->idose[j]){
           ind->ixds = j;
           break;
         }
       }
-    }
-    if (ind->ix[ind->idx] != ind->idose[ind->ixds]){
-      //524288
-      if (!(ind->err & rxErrSync2)){
-        ind->err += rxErrSync2;
+      if (ind->ix[ind->idx] != ind->idose[ind->ixds]){
+        for(int j = ind->ixds+1; j< ind->ndoses; j++){
+          if (ind->ix[ind->idx] == ind->idose[j]){
+            ind->ixds = j;
+            break;
+          }
+        }
+      }
+      if (ind->ix[ind->idx] != ind->idose[ind->ixds]){
+        //524288
+        if (!(ind->err & rxErrSync2)){
+          ind->err += rxErrSync2;
+        }
+        return 0;
       }
-      return 0;
     }
+    return 1;
   }
-  return 1;
-}
 
 
-extern t_F AMT;
+  extern t_F AMT;
 
 
-static inline double getAmt(rx_solving_options_ind *ind, int id, int cmt, double dose, double t, double *y) {
-  double ret = AMT(id, cmt, dose, t, y);
-  if (ISNA(ret)){
-    rx_solving_options *op = &op_global;
-    op->badSolve=1;
-    op->naTime = 1;
+  static inline double getAmt(rx_solving_options_ind *ind, int id, int cmt, double dose, double t, double *y) {
+    double ret = AMT(id, cmt, dose, t, y);
+    if (ISNA(ret)){
+      rx_solving_options *op = &op_global;
+      op->badSolve=1;
+      op->naTime = 1;
+    }
+    return ret;
   }
-  return ret;
-}
 
-static inline int isIgnoredDose(rx_solving_options_ind *ind, int ixds) {
-  for (int i = 0; i < ind->ignoredDosesN[0]; ++i) {
-    if (ind->idx < 0 ) {
-      return 0;
-    } else if (ind->idx >= 0 && ind->ignoredDoses[i] == ixds) {
-      return 1;
+  static inline int isIgnoredDose(rx_solving_options_ind *ind, int ixds) {
+    for (int i = 0; i < ind->ignoredDosesN[0]; ++i) {
+      if (ind->idx < 0 ) {
+        return 0;
+      } else if (ind->idx >= 0 && ind->ignoredDoses[i] == ixds) {
+        return 1;
+      }
     }
+    return 0;
   }
-  return 0;
-}
 
-static inline int pushIgnoredDose(int doseIdx, rx_solving_options_ind *ind) {
-  int re = 0;
-  for (int i = 0; i < ind->ignoredDosesN[0]; ++i) {
-    if (ind->ignoredDoses[i] == doseIdx) return 0;
-  }
-  if (ind->ignoredDosesN[0]+1 >= ind->ignoredDosesAllocN[0]) {
-    int *tmpI = (int*)realloc(ind->ignoredDoses, (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
-    if (tmpI == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 0;
+  static inline int pushIgnoredDose(int doseIdx, rx_solving_options_ind *ind) {
+    int re = 0;
+    for (int i = 0; i < ind->ignoredDosesN[0]; ++i) {
+      if (ind->ignoredDoses[i] == doseIdx) return 0;
     }
-    ind->ignoredDoses = tmpI;
-    ind->ignoredDosesAllocN[0] = (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE);
-    re = 1;
+    if (ind->ignoredDosesN[0]+1 >= ind->ignoredDosesAllocN[0]) {
+      int *tmpI = (int*)realloc(ind->ignoredDoses, (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
+      if (tmpI == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 0;
+      }
+      ind->ignoredDoses = tmpI;
+      ind->ignoredDosesAllocN[0] = (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE);
+      re = 1;
+    }
+    ind->ignoredDoses[ind->ignoredDosesN[0]] = doseIdx;
+    ind->ignoredDosesN[0] = ind->ignoredDosesN[0]+1;
+    return re;
   }
-  ind->ignoredDoses[ind->ignoredDosesN[0]] = doseIdx;
-  ind->ignoredDosesN[0] = ind->ignoredDosesN[0]+1;
-  return re;
-}
 
-static inline int pushPendingDose(int doseIdx, rx_solving_options_ind *ind) {
-  int re = 0;
-  if (ind->pendingDosesN[0]+1 >= ind->pendingDosesAllocN[0]) {
-    int *tmpI = (int*)realloc(ind->pendingDoses, (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
-    if (tmpI == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 0;
+  static inline int pushPendingDose(int doseIdx, rx_solving_options_ind *ind) {
+    int re = 0;
+    if (ind->pendingDosesN[0]+1 >= ind->pendingDosesAllocN[0]) {
+      int *tmpI = (int*)realloc(ind->pendingDoses, (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
+      if (tmpI == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 0;
+      }
+      ind->pendingDoses = tmpI;
+      ind->pendingDosesAllocN[0] = (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE);
+      re = 1;
     }
-    ind->pendingDoses = tmpI;
-    ind->pendingDosesAllocN[0] = (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE);
-    re = 1;
+    ind->pendingDoses[ind->pendingDosesN[0]] = doseIdx;
+    ind->pendingDosesN[0] = ind->pendingDosesN[0]+1;
+    return re;
   }
-  ind->pendingDoses[ind->pendingDosesN[0]] = doseIdx;
-  ind->pendingDosesN[0] = ind->pendingDosesN[0]+1;
-  return re;
-}
-
 
-static inline int pushDosingEvent(double time, double amt, int evid,
-                                   rx_solving_options_ind *ind) {
-  int re = 0;
-  if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
-    int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
-    if (tmpI == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 0;
-    }
-    ind->extraDoseTimeIdx = tmpI;
 
-    tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
-    if (tmpI == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 1;
-    }
-    ind->extraDoseEvid = tmpI;
+  static inline int pushDosingEvent(double time, double amt, int evid,
+                                    rx_solving_options_ind *ind) {
+    int re = 0;
+    if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
+      int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
+      if (tmpI == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 0;
+      }
+      ind->extraDoseTimeIdx = tmpI;
 
-    double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
-    if (tmpD == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 1;
-    }
-    ind->extraDoseTime = tmpD;
+      tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
+      if (tmpI == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 1;
+      }
+      ind->extraDoseEvid = tmpI;
 
-    tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
-    if (tmpD == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 1;
-    }
-    ind->extraDoseDose = tmpD;
+      double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
+      if (tmpD == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 1;
+      }
+      ind->extraDoseTime = tmpD;
 
-    ind->extraDoseAllocN[0] = (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE);
-    re = 1;
-  }
-  ind->extraDoseTimeIdx[ind->extraDoseN[0]] = ind->extraDoseN[0];
-  ind->extraDoseTime[ind->extraDoseN[0]] = time;
-  ind->extraDoseDose[ind->extraDoseN[0]] = amt;
-  ind->extraDoseEvid[ind->extraDoseN[0]] = evid;
-  pushPendingDose(-1-ind->extraDoseTimeIdx[ind->extraDoseN[0]], ind);
-  ind->extraDoseN[0] = ind->extraDoseN[0]+1;
-  ind->extraSorted = 0;
-  return re;
-}
+      tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
+      if (tmpD == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 1;
+      }
+      ind->extraDoseDose = tmpD;
 
-static inline int pushUniqueDosingEvent(double time, double amt, int evid,
-                                        rx_solving_options_ind *ind) {
-  int re = 0;
-  if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
-    int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
-    if (tmpI == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 0;
+      ind->extraDoseAllocN[0] = (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE);
+      re = 1;
     }
-    ind->extraDoseTimeIdx = tmpI;
+    ind->extraDoseTimeIdx[ind->extraDoseN[0]] = ind->extraDoseN[0];
+    ind->extraDoseTime[ind->extraDoseN[0]] = time;
+    ind->extraDoseDose[ind->extraDoseN[0]] = amt;
+    ind->extraDoseEvid[ind->extraDoseN[0]] = evid;
+    pushPendingDose(-1-ind->extraDoseTimeIdx[ind->extraDoseN[0]], ind);
+    ind->extraDoseN[0] = ind->extraDoseN[0]+1;
+    ind->extraSorted = 0;
+    return re;
+  }
 
-    tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
-    if (tmpI == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 1;
-    }
-    ind->extraDoseEvid = tmpI;
+  static inline int pushUniqueDosingEvent(double time, double amt, int evid,
+                                          rx_solving_options_ind *ind) {
+    int re = 0;
+    if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
+      int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
+      if (tmpI == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 0;
+      }
+      ind->extraDoseTimeIdx = tmpI;
 
-    double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
-    if (tmpD == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 1;
-    }
-    ind->extraDoseTime = tmpD;
+      tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
+      if (tmpI == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 1;
+      }
+      ind->extraDoseEvid = tmpI;
 
-    tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
-    if (tmpD == NULL) {
-      rx_solving_options *op = &op_global;
-      op->badSolve = 1;
-      return 1;
-    }
-    ind->extraDoseDose = tmpD;
+      double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
+      if (tmpD == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 1;
+      }
+      ind->extraDoseTime = tmpD;
 
-    ind->extraDoseAllocN[0] = (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE);
-    re = 1;
-  }
-  for (int i = 0; i < ind->extraDoseN[0]; ++i) {
-    if (ind->extraDoseTime[i] == time &&
-        ind->extraDoseDose[i] == amt &&
-        ind->extraDoseEvid[i] == evid) {
-      // found return early
-      return re;
-    }
-  }
-  ind->extraDoseTimeIdx[ind->extraDoseN[0]] = ind->extraDoseN[0];
-  ind->extraDoseTime[ind->extraDoseN[0]] = time;
-  ind->extraDoseDose[ind->extraDoseN[0]] = amt;
-  ind->extraDoseEvid[ind->extraDoseN[0]] = evid;
-  pushPendingDose(-1-ind->extraDoseTimeIdx[ind->extraDoseN[0]], ind);
-  ind->extraDoseN[0] = ind->extraDoseN[0]+1;
-  ind->extraSorted = 0;
-  return re;
-}
+      tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
+      if (tmpD == NULL) {
+        rx_solving_options *op = &op_global;
+        op->badSolve = 1;
+        return 1;
+      }
+      ind->extraDoseDose = tmpD;
 
-static inline int handle_evid(int evid, double *yp, double xout, rx_solving_options_ind *ind) {
-  rx_solving_options *op = &op_global;
-  int id  = ind->solveid;
-  if (isObs(evid)) return 0;
-  if (isIgnoredDose(ind, ind->ixds)) {
-    // REprintf("ignored evid %d dose at time %f is value %f (ind->ixds: %d: ind->idx: %d)\n",
-    //          evid, xout, getDoseIndex(ind, ind->idx), ind->ixds, ind->idx);
-    ind->ixds++;
-    ind->solved = ind->idx;
-    return 0;
-  }//  else if (!ind->doSS) {
-  //   REprintf("handle evid %d dose at time %f is value %f (ind->ixds: %d; ind->idx: %d)\n",
-  //            evid, xout, getDoseIndex(ind, ind->idx), ind->ixds, ind->idx);
-  // }
-  int cmt, foundBad, j;
-  double tmp;
-  getWh(evid, &(ind->wh), &(ind->cmt), &(ind->wh100), &(ind->whI), &(ind->wh0));
-  handleTlastInline(&xout, ind);
-  if (ind->wh0 == EVID0_SSINF) {
-    ind->ixds++;
-    ind->solved = ind->idx;
-    return 1;
-  }
-  /* wh100 = ind->wh100; */
-  cmt = ind->cmt;
-  if (cmt<0) {
-    if (!(ind->err & rxErrNegCmt)) {
-      ind->err += rxErrNegCmt;
-      /* Rprintf("Supplied an invalid EVID (EVID=%d; cmt %d)", evid, cmt); */
+      ind->extraDoseAllocN[0] = (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE);
+      re = 1;
     }
-    return 0;
-  }
-  if (cmt >= op->neq + op->extraCmt) {
-    foundBad = 0;
-    for (j = 0; j < ind->nBadDose; j++) {
-      if (ind->BadDose[j] == cmt+1) {
-        foundBad=1;
-        break;
+    for (int i = 0; i < ind->extraDoseN[0]; ++i) {
+      if (ind->extraDoseTime[i] == time &&
+          ind->extraDoseDose[i] == amt &&
+          ind->extraDoseEvid[i] == evid) {
+        // found return early
+        return re;
       }
     }
-    if (!foundBad) {
-      ind->BadDose[ind->nBadDose]=cmt+1;
-      ind->nBadDose++;
-    }
-  } else {
-    //if (syncIdx(ind) == 0) return 0;
-    if (ind->wh0 == EVID0_OFF) {
-      yp[cmt]=op_global.inits[cmt];
-      ind->InfusionRate[cmt] = 0;
-      ind->cacheME=0;
-      ind->on[cmt] = 0;
+    ind->extraDoseTimeIdx[ind->extraDoseN[0]] = ind->extraDoseN[0];
+    ind->extraDoseTime[ind->extraDoseN[0]] = time;
+    ind->extraDoseDose[ind->extraDoseN[0]] = amt;
+    ind->extraDoseEvid[ind->extraDoseN[0]] = evid;
+    pushPendingDose(-1-ind->extraDoseTimeIdx[ind->extraDoseN[0]], ind);
+    ind->extraDoseN[0] = ind->extraDoseN[0]+1;
+    ind->extraSorted = 0;
+    return re;
+  }
+
+  static inline int handle_evid(int evid, double *yp, double xout, rx_solving_options_ind *ind) {
+    rx_solving_options *op = &op_global;
+    int id  = ind->solveid;
+    if (isObs(evid)) return 0;
+    if (isIgnoredDose(ind, ind->ixds)) {
+      // REprintf("ignored evid %d dose at time %f is value %f (ind->ixds: %d: ind->idx: %d)\n",
+      //          evid, xout, getDoseIndex(ind, ind->idx), ind->ixds, ind->idx);
+      ind->ixds++;
+      ind->solved = ind->idx;
+      return 0;
+    }//  else if (!ind->doSS) {
+    //   REprintf("handle evid %d dose at time %f is value %f (ind->ixds: %d; ind->idx: %d)\n",
+    //            evid, xout, getDoseIndex(ind, ind->idx), ind->ixds, ind->idx);
+    // }
+    int cmt, foundBad, j;
+    double tmp;
+    getWh(evid, &(ind->wh), &(ind->cmt), &(ind->wh100), &(ind->whI), &(ind->wh0));
+    handleTlastInline(&xout, ind);
+    if (ind->wh0 == EVID0_SSINF) {
+      ind->ixds++;
+      ind->solved = ind->idx;
       return 1;
     }
-    if (!ind->doSS && (ind->wh0 == EVID0_SS2 || ind->wh0 == EVID0_SS20)) {
-      // Save for adding at the end; Only for ODE systems
-      memcpy(ind->solveSave, yp, (op->neq+op->extraCmt)*sizeof(double));
+    /* wh100 = ind->wh100; */
+    cmt = ind->cmt;
+    if (cmt<0) {
+      if (!(ind->err & rxErrNegCmt)) {
+        ind->err += rxErrNegCmt;
+        /* Rprintf("Supplied an invalid EVID (EVID=%d; cmt %d)", evid, cmt); */
+      }
+      return 0;
     }
-    switch(ind->whI) {
-    case EVIDF_MODEL_RATE_ON: // modeled rate.
-    case EVIDF_MODEL_DUR_ON: // modeled duration.
-      // Rate already calculated and saved in the next dose record
-      if (ind->wh0 != EVID0_SS0 &&
-          ind->wh0 != EVID0_SS20) {
-        ind->on[cmt] = 1;
-        ind->cacheME = 0;
-        ind->InfusionRate[cmt] -= getDoseIndexPlus1(ind, ind->idx);
-        // if (ind->wh0 != EVID0_SS2 &&
-        //     ind->wh0 != EVID0_SS) {
-        //   int infEixds = ind->ixds;
-        //   if (infEixds > 0) {
-        //     pushPendingDose(infEixds+1, ind);
-        //   } else {
-        //     pushPendingDose(infEixds-1, ind);
-        //   }
-        // }
+    if (cmt >= op->neq + op->extraCmt) {
+      foundBad = 0;
+      for (j = 0; j < ind->nBadDose; j++) {
+        if (ind->BadDose[j] == cmt+1) {
+          foundBad=1;
+          break;
+        }
+      }
+      if (!foundBad) {
+        ind->BadDose[ind->nBadDose]=cmt+1;
+        ind->nBadDose++;
+      }
+    } else {
+      //if (syncIdx(ind) == 0) return 0;
+      if (ind->wh0 == EVID0_OFF) {
+        yp[cmt]=op_global.inits[cmt];
+        ind->InfusionRate[cmt] = 0;
+        ind->cacheME=0;
+        ind->on[cmt] = 0;
+        return 1;
+      }
+      if (!ind->doSS && (ind->wh0 == EVID0_SS2 || ind->wh0 == EVID0_SS20)) {
+        // Save for adding at the end; Only for ODE systems
+        memcpy(ind->solveSave, yp, (op->neq+op->extraCmt)*sizeof(double));
+      }
+      switch(ind->whI) {
+      case EVIDF_MODEL_RATE_ON: // modeled rate.
+      case EVIDF_MODEL_DUR_ON: // modeled duration.
+        // Rate already calculated and saved in the next dose record
+        if (ind->wh0 != EVID0_SS0 &&
+            ind->wh0 != EVID0_SS20) {
+          ind->on[cmt] = 1;
+          ind->cacheME = 0;
+          ind->InfusionRate[cmt] -= getDoseIndexPlus1(ind, ind->idx);
+          // if (ind->wh0 != EVID0_SS2 &&
+          //     ind->wh0 != EVID0_SS) {
+          //   int infEixds = ind->ixds;
+          //   if (infEixds > 0) {
+          //     pushPendingDose(infEixds+1, ind);
+          //   } else {
+          //     pushPendingDose(infEixds-1, ind);
+          //   }
+          // }
+          if (ind->wh0 == EVID0_SS2 &&
+              getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp) !=
+              getDoseIndex(ind, ind->idx)) {
+            if (!(ind->err & rxErrModeledFss2)){
+              ind->err += rxErrModeledFss2;
+            }
+            return 0;
+          }
+        }
+        break;
+      case EVIDF_MODEL_RATE_OFF: // End modeled rate
+      case EVIDF_MODEL_DUR_OFF: // end modeled duration
+        // In this case re-sort is not going to be assessed
+        // If cmt is off, don't remove rate....
+        // Probably should throw an error if the infusion rate is on still.
+        // ind->curDose and ind->curDoseS[cmt] are handled when the modeled item is turned on.
+        ind->InfusionRate[cmt] += getDoseIndex(ind, ind->idx);
+        ind->cacheME=0;
         if (ind->wh0 == EVID0_SS2 &&
             getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp) !=
             getDoseIndex(ind, ind->idx)) {
-          if (!(ind->err & rxErrModeledFss2)){
-            ind->err += rxErrModeledFss2;
+          if (!(ind->err & 2097152)){
+            ind->err += 2097152;
           }
           return 0;
         }
-      }
-      break;
-    case EVIDF_MODEL_RATE_OFF: // End modeled rate
-    case EVIDF_MODEL_DUR_OFF: // end modeled duration
-      // In this case re-sort is not going to be assessed
-      // If cmt is off, don't remove rate....
-      // Probably should throw an error if the infusion rate is on still.
-      // ind->curDose and ind->curDoseS[cmt] are handled when the modeled item is turned on.
-      ind->InfusionRate[cmt] += getDoseIndex(ind, ind->idx);
-      ind->cacheME=0;
-      if (ind->wh0 == EVID0_SS2 &&
-          getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp) !=
-          getDoseIndex(ind, ind->idx)) {
-        if (!(ind->err & 2097152)){
-          ind->err += 2097152;
+        break;
+      case EVIDF_INF_DUR:
+        // In this case bio-availability changes the rate, but the
+        // duration remains constant.  rate = amt/dur
+        ind->on[cmt] = 1;
+        tmp = getDoseIndex(ind, ind->idx);
+        if (tmp > 0) {
+          ind->curDose = tmp;
+          ind->curDoseS[cmt] = ind->curDose;
+          // if (ind->wh0 != EVID0_SS2 &&
+          //     ind->wh0 != EVID0_SS) {
+          //   int infEixds;
+          //   handleInfusionGetEndOfInfusionIndex(ind->ixds, &infEixds, &rx_global, op, ind);
+          //   pushPendingDose(infEixds, ind);
+          // }
         }
-        return 0;
-      }
-      break;
-    case EVIDF_INF_DUR:
-      // In this case bio-availability changes the rate, but the
-      // duration remains constant.  rate = amt/dur
-      ind->on[cmt] = 1;
-      tmp = getDoseIndex(ind, ind->idx);
-      if (tmp > 0) {
-        ind->curDose = tmp;
-        ind->curDoseS[cmt] = ind->curDose;
-        // if (ind->wh0 != EVID0_SS2 &&
-        //     ind->wh0 != EVID0_SS) {
-        //   int infEixds;
-        //   handleInfusionGetEndOfInfusionIndex(ind->ixds, &infEixds, &rx_global, op, ind);
-        //   pushPendingDose(infEixds, ind);
-        // }
-      }
-      tmp = getAmt(ind, id, cmt, tmp, xout, yp);
-      ind->InfusionRate[cmt] += tmp;
-      ind->cacheME=0;
-      if (ind->wh0 == EVID0_SS2 && tmp != getDoseIndex(ind, ind->idx)) {
-        if (!(ind->err & 4194304)){
-          ind->err += 4194304;
+        tmp = getAmt(ind, id, cmt, tmp, xout, yp);
+        ind->InfusionRate[cmt] += tmp;
+        ind->cacheME=0;
+        if (ind->wh0 == EVID0_SS2 && tmp != getDoseIndex(ind, ind->idx)) {
+          if (!(ind->err & 4194304)){
+            ind->err += 4194304;
+          }
+          return 0;
         }
-        return 0;
-      }
-      break;
-    case EVIDF_INF_RATE:
-      // In this case bio-availability changes the duration, but the
-      // rate remains constant.  rate = amt/dur
-      ind->on[cmt] = 1;
-      tmp = getDoseIndex(ind, ind->idx);
-      if (tmp > 0) {
-        ind->curDose = tmp;
-        ind->curDoseS[cmt] = ind->curDose;
-        // if (ind->wh0 != EVID0_SS2 &&
-        //     ind->wh0 != EVID0_SS) {
-        //   int infEixds;
-        //   handleInfusionGetEndOfInfusionIndex(ind->ixds, &infEixds, &rx_global, op, ind);
-        //   pushPendingDose(infEixds, ind);
+        break;
+      case EVIDF_INF_RATE:
+        // In this case bio-availability changes the duration, but the
+        // rate remains constant.  rate = amt/dur
+        ind->on[cmt] = 1;
+        tmp = getDoseIndex(ind, ind->idx);
+        if (tmp > 0) {
+          ind->curDose = tmp;
+          ind->curDoseS[cmt] = ind->curDose;
+          // if (ind->wh0 != EVID0_SS2 &&
+          //     ind->wh0 != EVID0_SS) {
+          //   int infEixds;
+          //   handleInfusionGetEndOfInfusionIndex(ind->ixds, &infEixds, &rx_global, op, ind);
+          //   pushPendingDose(infEixds, ind);
+          // }
+        }
+        // if (!ind->doSS) {
+        //   REprintf("infusion dose at %f is %f ind->ixds: %d\n", xout, tmp, ind->ixds);
         // }
-      }
-      // if (!ind->doSS) {
-      //   REprintf("infusion dose at %f is %f ind->ixds: %d\n", xout, tmp, ind->ixds);
-      // }
-      ind->InfusionRate[cmt] += tmp;
-      ind->cacheME=0;
-      if (ind->wh0 == EVID0_SS2 && getDoseIndex(ind, ind->idx) > 0 &&
-          getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp) !=
-          getDoseIndex(ind, ind->idx)) {
-        if (!(ind->err & 4194304)){
-          ind->err += 4194304;
+        ind->InfusionRate[cmt] += tmp;
+        ind->cacheME=0;
+        if (ind->wh0 == EVID0_SS2 && getDoseIndex(ind, ind->idx) > 0 &&
+            getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp) !=
+            getDoseIndex(ind, ind->idx)) {
+          if (!(ind->err & 4194304)){
+            ind->err += 4194304;
+          }
+        }
+        break;
+      case EVIDF_REPLACE: // replace
+        ind->on[cmt] = 1;
+        yp[cmt] = getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp);     //dosing before obs
+        break;
+      case EVIDF_MULT: //multiply
+        ind->on[cmt] = 1;
+        yp[cmt] *= getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp);     //dosing before obs
+        break;
+      case EVIDF_NORMAL:
+        ind->on[cmt] = 1;
+        if (ind->wh0 != EVID0_PHANTOM) {
+          yp[cmt] += getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp);     //dosing before obs
         }
       }
-      break;
-    case EVIDF_REPLACE: // replace
-      ind->on[cmt] = 1;
-      yp[cmt] = getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp);     //dosing before obs
-      break;
-    case EVIDF_MULT: //multiply
-      ind->on[cmt] = 1;
-      yp[cmt] *= getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp);     //dosing before obs
-      break;
-    case EVIDF_NORMAL:
-      ind->on[cmt] = 1;
-      if (ind->wh0 != EVID0_PHANTOM) {
-        yp[cmt] += getAmt(ind, id, cmt, getDoseIndex(ind, ind->idx), xout, yp);     //dosing before obs
-      }
-		}
-		ind->ixds++;
-		ind->solved = ind->idx;
-    return 1;
-	}
-  return 0;
-}
+      ind->ixds++;
+      ind->solved = ind->idx;
+      return 1;
+    }
+    return 0;
+  }
 
-static inline int handleEvid1(int *i, rx_solve *rx, int *neq, double *yp, double *xout) {
-  rx_solving_options_ind *ind = &(rx->subjects[neq[1]]);
-  ind->idx = *i;
-  if (!isObs(getEvid(ind, ind->ix[ind->idx]))) {
-    syncIdx(ind);
+  static inline int handleEvid1(int *i, rx_solve *rx, int *neq, double *yp, double *xout) {
+    rx_solving_options_ind *ind = &(rx->subjects[neq[1]]);
+    ind->idx = *i;
+    if (!isObs(getEvid(ind, ind->ix[ind->idx]))) {
+      syncIdx(ind);
+    }
+    return handle_evid(getEvid(ind, ind->ix[ind->idx]), yp, *xout, ind);
   }
-  return handle_evid(getEvid(ind, ind->ix[ind->idx]), yp, *xout, ind);
-}
 
-// time   amt rate          ii  addl evid            ss
-static inline int getEvidFlag(int cmt, double amt, double rate, double ii, int evid, double ss) {
-	// #define  30
-	if (evid == 7) {
-		if (cmt > 0) {
-			return EVID0_PHANTOM;
-		} else {
-			return -1; // bad phantom
-		}
-	}
-	if (ss == 1.0) {
-		if (ii > 0.0) {
-			if (cmt > 0) {
-				return EVID0_SS;
-			} else {
-				return -2; // bad steady state 1
-			}
-		}
-		if (ii == 0.0 && amt == 0.0) {
-			if (cmt > 0) {
-				return EVID0_SSINF;
-			} else {
-				return -3; // bad infinite steady state infusion
-			}
-		}
-	} else if (ss == 2.0 && ii > 0.0) {
-		if (cmt > 0) {
-			return EVID0_SS2;
-		} else {
-			return -4;
-		}
-	}
-	if (cmt < 0) {
-		// turn off the compartment
-		return EVID0_OFF;
-	}
-	if (ss == 0.0 || ISNA(ss)) return EVID0_REGULAR;
-	return EVID0_REGULAR;
-}
+  // time   amt rate          ii  addl evid            ss
+  static inline int getEvidFlag(int cmt, double amt, double rate, double ii, int evid, double ss) {
+    // #define  30
+    if (evid == 7) {
+      if (cmt > 0) {
+        return EVID0_PHANTOM;
+      } else {
+        return -1; // bad phantom
+      }
+    }
+    if (ss == 1.0) {
+      if (ii > 0.0) {
+        if (cmt > 0) {
+          return EVID0_SS;
+        } else {
+          return -2; // bad steady state 1
+        }
+      }
+      if (ii == 0.0 && amt == 0.0) {
+        if (cmt > 0) {
+          return EVID0_SSINF;
+        } else {
+          return -3; // bad infinite steady state infusion
+        }
+      }
+    } else if (ss == 2.0 && ii > 0.0) {
+      if (cmt > 0) {
+        return EVID0_SS2;
+      } else {
+        return -4;
+      }
+    }
+    if (cmt < 0) {
+      // turn off the compartment
+      return EVID0_OFF;
+    }
+    if (ss == 0.0 || ISNA(ss)) return EVID0_REGULAR;
+    return EVID0_REGULAR;
+  }
 
 
-static inline int getEvidRateI(int cmt, double amt, double rate, double dur, double ii, int evid, double ss) {
-	if (evid == 1) {
-		if (dur == 0.0) {
-			if (rate == -1.0) {
-				return EVIDF_MODEL_RATE_ON;
-				// #define EVIDF_MODEL_RATE_OFF 7
-			} else if  (rate == -2.0) {
-				return EVIDF_MODEL_DUR_ON;
-			} else if (rate > 0.0) {
-				return EVIDF_INF_RATE;
-				// #define EVIDF_MODEL_DUR_OFF  6
-			}
-		} else if (rate == 0.0) {
-			if (dur > 0.0 ) {
-				return EVIDF_INF_DUR;
-			}
-		}
-	} else if (evid == 5) {
-		// replace
-		return EVIDF_REPLACE;
-	} else if (evid == 6) {
-		return EVIDF_MULT;
-	}
-	return EVIDF_NORMAL;
-}
+  static inline int getEvidRateI(int cmt, double amt, double rate, double dur, double ii, int evid, double ss) {
+    if (evid == 1) {
+      if (dur == 0.0) {
+        if (rate == -1.0) {
+          return EVIDF_MODEL_RATE_ON;
+          // #define EVIDF_MODEL_RATE_OFF 7
+        } else if  (rate == -2.0) {
+          return EVIDF_MODEL_DUR_ON;
+        } else if (rate > 0.0) {
+          return EVIDF_INF_RATE;
+          // #define EVIDF_MODEL_DUR_OFF  6
+        }
+      } else if (rate == 0.0) {
+        if (dur > 0.0 ) {
+          return EVIDF_INF_DUR;
+        }
+      }
+    } else if (evid == 5) {
+      // replace
+      return EVIDF_REPLACE;
+    } else if (evid == 6) {
+      return EVIDF_MULT;
+    }
+    return EVIDF_NORMAL;
+  }
 
-static inline int getEvidClassic(int cmt, double amt, double rate, double dur, double ii, int evid, double ss) {
-	if (isObs(evid)) return evid;
-	int cmtP = cmt;
-	int cmt100, cmt99, rateI, flg;
-	if (cmtP < 0) cmtP = -cmtP;
-	if (cmtP <= 99){
-		cmt100=0;
-		cmt99=cmtP;
-	} else {
-		cmt100=cmtP/100;
-		cmt99=cmtP-cmt100*100;
-	}
-	rateI = getEvidRateI(cmt, amt, rate, dur, ii, evid, ss);
-	flg = getEvidFlag(cmt, amt, rate, ii, evid, ss);
-	return cmt100*100000+rateI*10000+cmt99*100+flg;
-}
+  static inline int getEvidClassic(int cmt, double amt, double rate, double dur, double ii, int evid, double ss) {
+    if (isObs(evid)) return evid;
+    int cmtP = cmt;
+    int cmt100, cmt99, rateI, flg;
+    if (cmtP < 0) cmtP = -cmtP;
+    if (cmtP <= 99){
+      cmt100=0;
+      cmt99=cmtP;
+    } else {
+      cmt100=cmtP/100;
+      cmt99=cmtP-cmt100*100;
+    }
+    rateI = getEvidRateI(cmt, amt, rate, dur, ii, evid, ss);
+    flg = getEvidFlag(cmt, amt, rate, ii, evid, ss);
+    return cmt100*100000+rateI*10000+cmt99*100+flg;
+  }
 
 #if defined(__cplusplus)
 }