diff --git a/NAMESPACE b/NAMESPACE
index b4f615d18..f25be5128 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -286,6 +286,8 @@ export(.rxWithSink)
 export(.rxWithSinkBoth)
 export(.rxWithWd)
 export(.s3register)
+export(.solComp2)
+export(.solComp3)
 export(.symengineFs)
 export(.toClassicEvid)
 export(.udfCallFunArg)
diff --git a/R/linCmt.R b/R/linCmt.R
new file mode 100644
index 000000000..0a9c458fa
--- /dev/null
+++ b/R/linCmt.R
@@ -0,0 +1,56 @@
+#' Calculate the lambdas and coefficients of the two compartment model
+#'
+#' @param k10 elimination rate
+#' @param k12 rate from central to peripheral compartment
+#' @param k21 rate from peripheral to central compartment
+#' @param cpp boolean; call the c++ interface used for stan gradients
+#' @return List with `L` vector and matrices `C1` and `C2`
+#' @export
+#' @keywords internal
+#' @author Matthew L. Fidler based on `wnl` package/paper, implemented
+#'   in C/C++
+#' @examples
+#' .solComp2(k10=0.1, k12=3, k21=1)
+.solComp2 <- function(k10, k12, k21, cpp=FALSE) {
+  checkmate::assertNumeric(k10, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k12, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k21, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertLogical(cpp, len=1, any.missing = FALSE)
+  if (cpp) {
+    .ret <- .Call(`_rxode2_solComp2cpp`, k10, k12, k21)
+  } else {
+    .ret <- .Call(`_rxode2_solComp2`, k10, k12, k21)
+  }
+  if (is.null(.ret)) {
+    stop("roots must be distinct real values", call.=FALSE)
+  }
+  .ret
+}
+#' Calculate the lambdas and coefficients of the three compartment model
+#'
+#' @inheritParams .solComp2
+#' @param k13 rate from central to peripheral compartment #2
+#' @param k31 rate from peripheral compartment #2 to central
+#' @return List with `L` vector and matrices `C1`, `C2` and `C3`
+#' @export
+#' @author Matthew L. Fidler
+#' @keywords internal
+#' @examples
+#' .solComp3(k10=0.1, k12=3, k21=1, k13=2, k31=0.5)
+.solComp3 <- function(k10, k12, k21, k13, k31, cpp=FALSE) {
+  checkmate::assertNumeric(k10, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k12, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k21, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k13, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k31, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertLogical(cpp, len=1, any.missing = FALSE)
+  if (cpp) {
+    .ret <- .Call(`_rxode2_solComp3cpp`, k10, k12, k21, k13, k31)
+  } else {
+    .ret <- .Call(`_rxode2_solComp3`, k10, k12, k21, k13, k31)
+  }
+  if (is.null(.ret)) {
+    stop("roots must be distinct real values", call.=FALSE)
+  }
+  .ret
+}
diff --git a/inst/include/rxode2.h b/inst/include/rxode2.h
index 2ac86c68d..634224a8f 100644
--- a/inst/include/rxode2.h
+++ b/inst/include/rxode2.h
@@ -30,11 +30,14 @@
 #define getSolve(idx) ind->solve + (op->neq + op->nlin)*(idx)
 
 #ifdef _isrxode2_
-
+#define min2( a , b )  ( (a) < (b) ? (a) : (b) )
 #define max2( a , b )  ( (a) > (b) ? (a) : (b) )
 #define isSameTime(xout, xp) (fabs((xout)-(xp))  <= DBL_EPSILON*max2(fabs(xout),fabs(xp)))
+
 // use ~dop853 definition of same time
 #define isSameTimeDop(xout, xp) (0.1 * fabs((xout)-(xp)) <= fabs(xout) * 2.3E-16)
+#define isSameTimeOp(xout, xp) (op->stiff == 0 ? isSameTimeDop(xout, xp) : isSameTime(xout, xp))
+
 
 #else
 
diff --git a/man/dot-solComp2.Rd b/man/dot-solComp2.Rd
new file mode 100644
index 000000000..379f0c8d2
--- /dev/null
+++ b/man/dot-solComp2.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/linCmt.R
+\name{.solComp2}
+\alias{.solComp2}
+\title{Calculate the lambdas and coefficients of the two compartment model}
+\usage{
+.solComp2(k10, k12, k21)
+}
+\arguments{
+\item{k10}{elimination rate}
+
+\item{k12}{rate from central to peripheral compartment}
+
+\item{k21}{rate from peripheral to central compartment}
+}
+\value{
+List with \code{L} vector and matrices \code{C1} and \code{C2}
+}
+\description{
+Calculate the lambdas and coefficients of the two compartment model
+}
+\examples{
+.solComp2(k10=0.1, k12=3, k21=1)
+}
+\author{
+Matthew L. Fidler based on \code{wnl} package/paper, implemented
+in C/C++
+}
+\keyword{internal}
diff --git a/man/dot-solComp3.Rd b/man/dot-solComp3.Rd
new file mode 100644
index 000000000..6322f8d81
--- /dev/null
+++ b/man/dot-solComp3.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/linCmt.R
+\name{.solComp3}
+\alias{.solComp3}
+\title{Calculate the lambdas and coefficients of the three compartment model}
+\usage{
+.solComp3(k10, k12, k21, k13, k31)
+}
+\arguments{
+\item{k10}{elimination rate}
+
+\item{k12}{rate from central to peripheral compartment}
+
+\item{k21}{rate from peripheral to central compartment}
+
+\item{k13}{rate from central to peripheral compartment #2}
+
+\item{k31}{rate from peripheral compartment #2 to central}
+}
+\value{
+List with \code{L} vector and matrices \code{C1}, \code{C2} and \code{C3}
+}
+\description{
+Calculate the lambdas and coefficients of the three compartment model
+}
+\examples{
+.solComp3(k10=0.1, k12=3, k21=1, k13=2, k31=0.5)
+}
+\author{
+Matthew L. Fidler
+}
+\keyword{internal}
diff --git a/src/approx.c b/src/approx.c
index f4f13b66d..a3f3bfc11 100644
--- a/src/approx.c
+++ b/src/approx.c
@@ -13,8 +13,6 @@
 #define _as_zero(a) (fabs(a) < sqrt(DBL_EPSILON) ? 0.0 : a)
 #define _as_dbleps(a) (fabs(a) < sqrt(DBL_EPSILON) ? ((a) < 0 ? -sqrt(DBL_EPSILON)  : sqrt(DBL_EPSILON)) : a)
 
-#define isSameTimeOp(xout, xp) (op->stiff == 0 ? isSameTimeDop(xout, xp) : isSameTime(xout, xp))
-
 #ifdef ENABLE_NLS
 #include <libintl.h>
 #define _(String) dgettext ("rxode2", String)
diff --git a/src/comp.c b/src/comp.c
new file mode 100644
index 000000000..12868ca37
--- /dev/null
+++ b/src/comp.c
@@ -0,0 +1,71 @@
+#define USE_FC_LEN_T
+#define STRICT_R_HEADERS
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>   /* dj: import intptr_t */
+//#include "ode.h"
+#define R_NO_REMAP
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
+#include <Rmath.h>
+#include <R_ext/Lapack.h>
+#include <R_ext/BLAS.h>
+#ifndef FCONE
+# define FCONE
+#endif
+
+#include "../inst/include/rxode2.h"
+#include "comp.h"
+#include "parTrans.h"
+
+#ifdef ENABLE_NLS
+#include <libintl.h>
+#define _(String) dgettext ("rxode2", String)
+/* replace pkg as appropriate */
+#else
+#define _(String) (String)
+#endif
+
+int comp1c(int ncmt,
+           int trans,
+           double *yp, // prior y
+           double xp, // initial time point
+           double xout, // final time point
+           double *p,
+           double *rate) {
+  // no update for same time
+  if (isSameTime(xout, xp)) return 1;
+  lin_context_c_t lin;
+  lin.ncmt = ncmt;
+  lin.ka=lin.k10=
+    lin.k12=lin.k21=
+    lin.k13=lin.k31=lin.v=
+    lin.dt=xout-xp;
+  lin.rate = rate;
+  if (!parTrans(&trans, &p[0], &p[1], &p[2], &p[3], &p[4], &p[5],
+                &ncmt, &(lin.k10), &(lin.v), &(lin.k12),
+                &(lin.k21), &(lin.k13), &(lin.k31))){
+    return 0;
+  }
+  return comp1(yp, &lin);
+}
+
+// This is a R wrapper for comp1c
+SEXP _rxode2_comp1c(SEXP prior, SEXP xpS, SEXP xoutS,
+                    SEXP rateS, SEXP transS, SEXP parS) {
+  SEXP retS = PROTECT(Rf_allocVector(INTSXP, 1));
+  int *ret = INTEGER(retS);
+  ret[0] = 0;
+  int ncmt = INTEGER(VECTOR_ELT(transS, 1))[0];
+  int trans = INTEGER(VECTOR_ELT(transS, 2))[0];
+  comp1c(ncmt, trans, REAL(prior), // prior y
+         REAL(xpS)[0], // initial time point
+         REAL(xoutS)[0], // final time point
+         REAL(rateS),
+         REAL(rateS));
+  UNPROTECT(1);
+  return retS;
+}
diff --git a/src/comp.h b/src/comp.h
new file mode 100644
index 000000000..9ae8c2f31
--- /dev/null
+++ b/src/comp.h
@@ -0,0 +1,208 @@
+#ifndef __COMP_H__
+#define __COMP_H__
+
+#include "../inst/include/rxode2.h"
+#include "solComp.h"
+
+typedef struct lin_context_c_t {
+  double *rate; // rate (depot + central or depot only)
+  int ncmt;
+  int oral0; // does the system have an oral dose? (0:central or 1:depot+central)
+  // typical PK constants and v for volume of distribution
+  double ka;
+  double k10;
+  double k12;
+  double k21;
+  double k13;
+  double k31;
+  double v;
+  // dt -- delta t
+  double dt;
+} lin_context_c_t;
+
+static inline void printSqMat(double *in, int d) {
+  // debugging function to print a square matrix in R.
+  int pro=0;
+  SEXP mat = PROTECT(Rf_allocVector(REALSXP, d*d)); pro++;
+  SEXP dm  = PROTECT(Rf_allocVector(INTSXP, 2)); pro++;
+  int *dmi = INTEGER(dm);
+  dmi[0] =dmi[1] = d;
+  double *matd = REAL(mat);
+  for (int i = d*d; i--;)  {
+    matd[i] = in[i];
+  }
+  Rf_setAttrib(mat, R_DimSymbol, dm);
+  Rf_PrintValue(mat);
+  UNPROTECT(pro);
+}
+
+static inline void printVec(double *in, int d) {
+  int pro=0;
+  SEXP vec = PROTECT(Rf_allocVector(REALSXP, d)); pro++;
+  double *vecd =REAL(vec);
+  for (int i = d; i--;) {
+    vecd[i] = in[i];
+  }
+  Rf_PrintValue(vec);
+  UNPROTECT(pro);
+}
+
+// Handle single point solve
+static inline int comp1solve1(double *yp, lin_context_c_t *lin) {
+  double E  = exp(-(lin->k10)*lin->dt);
+  double Ea = E;
+  double pDepot = 0.0;
+  double rDepot = 0.0;
+  double R = lin->rate[lin->oral0];
+  // In the derivation used in https://doi.org/10.12793/tcp.2019.27.2.43
+  // the expression for the Ka contribution (Eq 9, Eq 13 and Eq 33) is given by
+  //
+  // Ka*Xg(0)*exp(-Ka*t)
+  //
+  // This is true as long as there is not an infusion in the depot
+  //
+  // When there is an infusion (Rg) in the depot the ka would be:
+  //
+  // Ka*[Kg(0)*exp(-Ka*t) + Rg/ka*(1-exp(-Ka*t))]
+  //
+  // as implied by equation #12 (which is the eq in the bracket)
+  //
+  // expanding this becomes:
+  //
+  // (Ka*Kg(0) - Rg)*exp(-Ka*t) + Rg
+  //
+  // Both Ka*Kg(0) and Ka*kg(0)-Rg in general are not dependent on
+  // time.  Also Rg is simply a time invariant constant
+  //
+  // Which means to get equations where infusions into a depot are supported
+  // You need to simply change 2 items:
+  //
+  // - in Eq 11, 32 and 41 you change Ka*Kg(0) to (Ka*Kg(0)- Rg)
+  // - in Eq 11, 32 and 41 you change R to (R+Rg)
+  //
+  // This was observed after solving a few systems manually
+  if (lin->oral0 == 1) {
+    Ea = exp(-(lin->ka)*lin->dt);
+    pDepot = yp[0];
+    rDepot = lin->rate[0];
+    R = rDepot + R;
+  }
+  yp[lin->oral0] = yp[lin->oral0]*E + R*(1.0-E)/(lin->k10);
+  if (isSameTime(lin->ka, lin->k10)) {
+    yp[lin->oral0] += (pDepot*(lin->k10)-rDepot)*lin->dt*E;
+  } else {
+    yp[lin->oral0] += (pDepot*(lin->ka)-rDepot)*(E-Ea)/((lin->ka)-(lin->k10));
+  }
+  if (lin->oral0) {
+    yp[0] = rDepot*(1.0-Ea)/(lin->ka) + pDepot*Ea;
+  }
+  return 1;
+}
+
+// prior solving information, will be updated with new information
+// (like lsoda and the like)
+static inline int comp1solve2(double *yp, lin_context_c_t *lin) {
+  double L[2], C1[4], C2[4], E[2], Ea[2], Xo[2], Rm[2];
+  double rDepot=0.0;
+  double R=lin->rate[lin->oral0];
+  //  double dT = lin->dt;
+  if (solComp2C(&(lin->k10), &(lin->k12), &(lin->k21), L, C1, C2) == 0) {
+    return 0;
+  }
+  Xo[0] = Xo[1] = 0.0;
+  E[0] = Ea[0] = exp(-L[0]*lin->dt);
+  E[1] = Ea[1] = exp(-L[1]*lin->dt);
+  const double one = 1.0;
+  const int ione = 1,
+    itwo = 2;
+  //Xo = Xo + pX[1 + j] * Co[, , j] %*% E # Bolus
+  F77_CALL(dgemm)("N", "N", &itwo, &ione, &itwo, &(yp[lin->oral0]), C1, &itwo, E,
+                  &itwo, &one, Xo, &itwo FCONE FCONE);
+  F77_CALL(dgemm)("N", "N", &itwo, &ione, &itwo, &(yp[lin->oral0+1]), C2, &itwo, E,
+                  &itwo, &one, Xo, &itwo FCONE FCONE);
+  if (lin->oral0 == 1 && yp[0] > 0.0) {
+    // Xo = Xo + Ka*pX[1]*(Co[, , 1] %*% ((E - Ea)/(Ka - L)))
+    rDepot = lin->rate[0];
+    R += rDepot;
+    double expa = exp(-(lin->ka)*lin->dt);
+    Ea[0] = (E[0]- expa)/((lin->ka)-L[0]);
+    Ea[1] = (E[1]- expa)/((lin->ka)-L[1]);
+    double cf = (lin->ka)*yp[0] - rDepot;
+    F77_CALL(dgemm)("N", "N", &itwo, &ione, &itwo, &cf, C1, &itwo, Ea, &itwo, &one, Xo, &itwo FCONE FCONE);
+    yp[0] = rDepot*(1.0-expa)/(lin->ka) + yp[0]*expa;
+  }
+  if (!isSameTime(R, 0.0)) {
+    // Xo = Xo + ((cR*Co[, , 1]) %*% ((1 - E)/L)) # Infusion
+    Rm[0] = (1.0 - E[0])/L[0];
+    Rm[1] = (1.0 - E[1])/L[1];
+    F77_CALL(dgemm)("N", "N", &itwo, &ione, &itwo, &(R),
+                    C1, &itwo, Rm, &itwo, &one, Xo, &itwo FCONE FCONE);
+  }
+  yp[lin->oral0] = Xo[0];
+  yp[lin->oral0+1] = Xo[1];
+  return 1;
+}
+
+static inline int comp1solve3(double *yp, lin_context_c_t *lin) {
+  double L[3], C1[9], C2[9], C3[9], E[3], Ea[3], Xo[3], Rm[3];
+  double R = lin->rate[lin->oral0];
+  double rDepot = 0.0;
+  if (solComp3C(&(lin->k10), &(lin->k12), &(lin->k21),
+                &(lin->k13), &(lin->k31), L, C1, C2, C3) == 0) {
+    return 0;
+  }
+  E[0] = Ea[0] = exp(-L[0]*lin->dt);
+  E[1] = Ea[1] = exp(-L[1]*lin->dt);
+  E[2] = Ea[2] = exp(-L[2]*lin->dt);
+  const double one = 1.0;
+  const int ione = 1, ithree=3;
+  //Xo = Xo + pX[1 + j] * Co[, , j] %*% E # Bolus
+  Xo[0]=Xo[1]=Xo[2]=0.0;
+  F77_CALL(dgemm)("N", "N", &ithree, &ione, &ithree, &(yp[lin->oral0]), C1, &ithree,
+                  E, &ithree, &one, Xo, &ithree FCONE FCONE);
+  F77_CALL(dgemm)("N", "N", &ithree, &ione, &ithree, &(yp[lin->oral0+1]), C2, &ithree,
+                  E, &ithree, &one, Xo, &ithree FCONE FCONE);
+  F77_CALL(dgemm)("N", "N", &ithree, &ione, &ithree, &(yp[lin->oral0+2]), C3, &ithree,
+                  E, &ithree, &one, Xo, &ithree FCONE FCONE);
+  if (lin->oral0 == 1 && yp[0] > 0.0) {
+    // Xo = Xo + Ka*pX[1]*(Co[, , 1] %*% ((E - Ea)/(Ka - L)))
+    rDepot=lin->rate[0];
+    R += rDepot;
+    double expa = exp(-(lin->ka)*lin->dt);
+    Ea[0] = (E[0]- expa)/((lin->ka) - L[0]);
+    Ea[1] = (E[1]- expa)/((lin->ka) - L[1]);
+    Ea[2] = (E[2]- expa)/((lin->ka) - L[2]);
+    double expa2 = (lin->ka)*yp[0] - rDepot;
+    F77_CALL(dgemm)("N", "N", &ithree, &ione, &ithree, &expa2, C1, &ithree,
+                    Ea, &ithree, &one, Xo, &ithree FCONE FCONE);
+    yp[0] = rDepot*(1.0-expa)/(lin->ka) + yp[0]*expa;
+  }
+  if (!isSameTime(R, 0.0)) {
+    // Xo = Xo + ((cR*Co[, , 1]) %*% ((1 - E)/L)) # Infusion
+    Rm[0] = (1.0 - E[0])/L[0];
+    Rm[1] = (1.0 - E[1])/L[1];
+    Rm[2] = (1.0 - E[2])/L[2];
+    F77_CALL(dgemm)("N", "N", &ithree, &ione, &ithree, &(R), C1, &ithree,
+                    Rm, &ithree, &one, Xo, &ithree FCONE FCONE);
+  }
+  yp[lin->oral0]   = Xo[0];
+  yp[lin->oral0+1] = Xo[1];
+  yp[lin->oral0/+2] = Xo[2];
+  return 1;
+}
+
+static inline int comp1(double *yp, lin_context_c_t *lin) {
+  switch (lin->ncmt) {
+  case 3:
+    return comp1solve3(yp, lin);
+    break;
+  case 2:
+    return comp1solve2(yp, lin);
+    break;
+  case 1:
+    return comp1solve1(yp, lin);
+    break;
+  }
+  return 0;
+}
+#endif
diff --git a/src/compB.cpp b/src/compB.cpp
new file mode 100644
index 000000000..e9bfa279b
--- /dev/null
+++ b/src/compB.cpp
@@ -0,0 +1,423 @@
+#define USE_FC_LEN_T
+#define STRICT_R_HEADERS
+#define R_NO_REMAP
+#define USE_FC_LEN_T
+#define STRICT_R_HEADERS
+#include <stan/math.hpp>
+#include <RcppEigen.h>
+#include "../inst/include/rxode2.h"
+#ifdef ENABLE_NLS
+#include <libintl.h>
+#define _(String) dgettext ("rxode2", String)
+/* replace pkg as appropriate */
+#else
+#define _(String) (String)
+#endif
+
+extern "C" void RSprintf(const char *format, ...);
+namespace stan {
+  namespace math {
+    using std::exp;
+    using stan::math::exp;
+    using std::sqrt;
+    using stan::math::sqrt;
+    using std::pow;
+    using stan::math::pow;
+    using std::acos;
+    using stan::math::acos;
+    using std::cos;
+    using stan::math::cos;
+
+    template<class T>
+    Eigen::Matrix<T, Eigen::Dynamic, 2>
+    parTransB(const Eigen::Matrix<T, Eigen::Dynamic, 1>& p,
+                  const int& ncmt,
+                  const int& trans) {
+      Eigen::Matrix<T, Eigen::Dynamic, 2> g(ncmt,3);
+      T btemp, ctemp, dtemp;
+#define p1    p[0]
+#define v1    p[1]
+#define p2    p[2]
+#define p3    p[3]
+#define p4    p[4]
+#define p5    p[5]
+#define v     g(0, 0)
+#define k     g(0, 1)
+
+#define k12   g(1, 0)
+#define k23   g(1, 0)
+
+#define k21   g(1, 1)
+#define k32   g(1, 1)
+
+#define k13   g(2, 0)
+#define k24   g(2, 0)
+
+#define k31   g(2, 1)
+#define k42   g(2, 1)
+      switch (ncmt) {
+      case 3: { // 3 compartment model
+        switch (trans){
+        case 1: // cl v q vp
+          k = p1/v1; // k = CL/V
+          v = v1;
+          k12 = p2/v1; // k12 = Q/V
+          k21 = p2/p3; // k21 = Q/Vp
+          k13 = p4/v1; // k31 = Q2/V
+          k31 = p4/p5; // k31 = Q2/Vp2
+          break;
+        case 2: // k=(*p1) v=(*v1) k12=(*p2) k21=(*p3) k13=(*p4) k31=(*p5)
+          k = p1;
+          v = v1;
+          k12 = p2;
+          k21 = p3;
+          k13 = p4;
+          k31 = p5;
+          break;
+        case 11:
+#undef beta
+#define A (1/v1)
+#define B (p3)
+#define C (p5)
+#define alpha (p1)
+#define beta (p2)
+#define gamma (p4)
+          v=1/(A+B+C);
+          btemp = -(alpha*C + alpha*B + gamma*A + gamma*B + beta*A + beta*C)*v;
+          ctemp = (alpha*beta*C + alpha*gamma*B + beta*gamma*A)*v;
+          dtemp = sqrt(btemp*btemp-4*ctemp);
+          k21 = 0.5*(-btemp+dtemp);
+          k31 = 0.5*(-btemp-dtemp);
+          k   = alpha*beta*gamma/k21/k31;
+          k12 = ((beta*gamma + alpha*beta + alpha*gamma) -
+                 k21*(alpha+beta+gamma) - k * k31 + k21*k21)/(k31 - k21);
+          k13 = alpha + beta + gamma - (k + k12 + k21 + k31);
+          break;
+        case 10:
+#undef A
+#define A v1
+          v=1/(A+B+C);
+          btemp = -(alpha*C + alpha*B + gamma*A + gamma*B + beta*A + beta*C)*v;
+          ctemp = (alpha*beta*C + alpha*gamma*B + beta*gamma*A)*v;
+          dtemp = sqrt(btemp*btemp-4*ctemp);
+          k21 = 0.5*(-btemp+dtemp);
+          k31 = 0.5*(-btemp-dtemp);
+          k   = alpha*beta*gamma/k21/k31;
+          k12 = ((beta*gamma + alpha*beta + alpha*gamma) -
+                 k21*(alpha+beta+gamma) - k * k31 + k21*k21)/(k31 - k21);
+          k13 = alpha + beta + gamma - (k + k12 + k21 + k31);
+#undef A
+#undef B
+#undef C
+#undef alpha
+#undef beta
+#undef gamma
+#define beta Rf_beta
+          break;
+        }
+      } break;
+      case 2:{ // 2 compartment model
+        switch (trans){
+        case 1: // cl=(*p1) v=(*v1) q=(*p2) vp=(*p3)
+          k = p1/v1; // k = CL/V
+          v = v1;
+          k12 = p2/v1; // k12 = Q/V
+          k21 = p2/p3; // k21 = Q/Vp
+          break;
+        case 2: // k=(*p1), (*v1)=v k12=(*p2) k21=(*p3)
+          k = p1;
+          v = v1;
+          k12 = p2;
+          k21 = p3;
+          break;
+        case 3: // cl=(*p1) v=(*v1) q=(*p2) vss=(*p3)
+          k = p1/v1; // k = CL/V
+          v = v1;
+          k12 = p2/v1; // k12 = Q/V
+          k21 = p2/(p3-v1); // k21 = Q/(Vss-V)
+          break;
+        case 4: // alpha=(*p1) beta=(*p2) k21=(*p3)
+          v = v1;
+          k21 = p3;
+          k = p1*p2/k21; // (*p1) = alpha (*p2) = beta
+          k12 = p1 + p2 - k21 - k;
+          break;
+        case 5: // alpha=(*p1) beta=(*p2) aob=(*p3)
+          v=v1;
+          k21 = (p3*p2+p1)/(p3+1);
+          k = (p1*p2)/k21;
+          k12 = p1+ p2 - k21 - k;
+          break;
+        case 11: // A2 V, alpha=(*p1), beta=(*p2), k21
+#undef beta
+#define A (1/v1)
+#define B (p3)
+#define alpha (p1)
+#define beta (p2)
+          v   = 1/(A+B);
+          k21 = (A*beta + B*alpha)*v;
+          k   = alpha*beta/k21;
+          k12 = alpha+beta-k21-k;
+          break;
+        case 10: // A=(*v1), alpha=(*p1), beta=(*p2), B=(*p3)
+          // Convert to A (right now A=(*v1) or A=1/(*v1))
+#undef A
+#define A (v1)
+          v   = 1/(A + B);
+          k21 = (A*beta + B*alpha)*v;
+          k   = alpha*beta/k21;
+          k12 = alpha + beta - k21 - k;
+#undef A
+#undef B
+#undef alpha
+#undef beta
+#define beta Rf_beta
+          break;
+        default:
+          RSprintf(_("invalid trans (2 cmt trans %d)\n"), trans);
+          return g;
+        }
+      } break;
+      case 1:{ // One compartment model
+        switch(trans){
+        case 1: // cl v
+          k = p1/v1; // k = CL/V
+          v = v1;
+          break;
+        case 2: // k V
+          k = p1;
+          v = v1;
+          break;
+        case 11: // alpha V
+          k = p1;
+          v = v1;
+          break;
+        case 10: // alpha A
+          k = p1;
+          v = 1/v1;
+          break;
+        default:
+          return g;
+        }
+      } break;
+      }
+#undef p1
+#undef v1
+#undef p2
+#undef p3
+#undef p4
+#undef p4
+
+#undef k
+#undef v
+#undef k12
+#undef k21
+#undef k13
+#undef k31
+      return g;
+    }
+  }
+
+#define v     g(0, 0)
+#define k20   g(0, 1)
+#define kel   g(0, 1)
+#define k23   g(1, 0)
+#define k32   g(1, 1)
+#define k24   g(2, 0)
+#define k42   g(2, 1)
+#define k10   g(0, 1)
+#define k12   g(1, 0)
+#define k21   g(1, 1)
+#define k13   g(2, 0)
+#define k31   g(2, 1)
+
+  template<class T>
+  int
+  solComp2Cpp(const  Eigen::Matrix<T, Eigen::Dynamic, 2> g,
+              Eigen::Matrix<T, 2, 1>& L,
+              Eigen::Matrix<T, 2, 2>& C1,
+              Eigen::Matrix<T, 2, 2>& C2) {
+
+    Eigen::Matrix<T, 2, 1> div;
+
+    T sum = (k10) + (k12) + (k21);
+    T disc= sqrt(sum*sum - 4.0* (k10)*(k21));
+    T tmp;
+    L(0, 0) = 0.5*(sum + disc);
+    L(1, 0) = 0.5*(sum - disc);
+    div(0, 0) = L(1, 0) - L(0, 0);
+    div(1, 0) = L(0, 0) - L(1, 0);
+    if (div(0, 0)*div(1, 0) == 0) return 0;
+    // c[0] = (0, 0)
+    // c[1] = (1, 0)
+    // c[2] = (0, 1)
+    // c[3] = (1, 1)
+
+    C1(0, 0) = (k21) - L(0, 0);
+    C1(0, 1) = (k21) - L(1, 0);
+    C2(0, 0) = C2(0, 1) = (k21);
+    C1(1, 0) = C1(1, 1) = (k12);
+    tmp = (k10) + (k12);
+    C2(1, 0) = tmp - L(0, 0);
+    C2(1, 1) = tmp - L(1, 0);
+    C1(0, 0) = C1(0, 0)/div(0, 0);
+    C1(1, 0) = C1(1, 0)/div(0, 0);
+    C2(0, 0) = C2(0, 0)/div(0, 0);
+    C2(1, 0) = C2(1, 0)/div(0, 0);
+    C1(0, 1) = C1(0, 1)/div(1, 0);
+    C1(1, 1) = C1(1, 1)/div(1, 0);
+    C2(0, 1) = C2(0, 1)/div(1, 0);
+    C2(1, 1) = C2(1, 1)/div(1, 0);
+    return 1;
+  }
+
+  template<class T>
+  int
+  solComp3Cpp(const  Eigen::Matrix<T, Eigen::Dynamic, 2> g,
+              Eigen::Matrix<T, 3, 1>& L,
+              Eigen::Matrix<T, 3, 3>& C1,
+              Eigen::Matrix<T, 3, 3>& C2,
+              Eigen::Matrix<T, 3, 3>& C3) {
+
+    T A1 = k10 + k12 + k13 + k21 + k31;
+    T A2 = k10*k21 + k10*k31 + k12*k31 +
+      k13*k21 + k21*k31;
+    T A3 = k21*k31*k10;
+    T Q  = (A1*A1 - 3.0*A2)/9.0;
+    T RQ = 2.0*sqrt(Q);
+    T R  = (2.0*A1*A1*A1 - 9.0*A1*A2 + 27.0*A3)/54.0;
+    T M  = Q*Q*Q - R*R;
+    if (M < 0) return 0;//stop("Error: Not real roots.")
+    T Th = acos(8.0*R/(RQ*RQ*RQ));
+    L(0, 0) = RQ*cos(Th/3.0) + A1/3.0;
+    L(1, 0) = RQ*cos((Th + M_2PI)/3.0) + A1/3.0;
+    L(2, 0) = RQ*cos((Th + 2*M_2PI)/3.0) + A1/3.0;
+
+    Eigen::Matrix<T, 3, 1> D;
+    D(0, 0) = (L(1, 0) - L(0, 0))*(L(2, 0) - L(0, 0));
+    D(1, 0) = (L(0, 0) - L(1, 0))*(L(2, 0) - L(1, 0));
+    D(2, 0) = (L(0, 0) - L(2, 0))*(L(1, 0) - L(2, 0));
+    if (D(0, 0)*D(1, 0)*D(2, 0) == 0.0) return 0;
+
+    C1(0, 0) = (k21 - L(0, 0))*(k31 - L(0, 0));
+    C1(0, 1) = (k21 - L(1, 0))*(k31 - L(1, 0));
+    C1(0, 2) = (k21 - L(2, 0))*(k31 - L(2, 0));
+
+    C2(0, 0) = (k21)*(k31 - L(0, 0));
+    C2(0, 1) = (k21)*(k31 - L(1, 0));
+    C2(0, 2) = (k21)*(k31 - L(2, 0));
+
+    C3(0, 0) = (k31)*(k21 - L(0, 0));
+    C3(0, 1) = (k31)*(k21 - L(1, 0));
+    C3(0, 2) = (k31)*(k21 - L(2, 0));
+
+    C1(1, 0) = (k12)*(k31 - L(0, 0));
+    C1(1, 1) = (k12)*(k31 - L(1, 0));
+    C1(1, 2) = (k12)*(k31 - L(2, 0));
+
+    C2(1, 0) = (k10 + k12 + k13 - L(0, 0))*(k31 - L(0, 0)) - (k31)*(k13);
+    C2(1, 1) = (k10 + k12 + k13 - L(1, 0))*(k31 - L(1, 0)) - (k31)*(k13);
+    C2(1, 2) = (k10 + k12 + k13 - L(2, 0))*(k31 - L(2, 0)) - (k31)*(k13);
+
+    C3(1, 0) = C3(1, 1) = C3(1, 2) = (k12)*(k31);
+
+    C1(2, 0) = (k13)*(k21 - L(0, 0));
+    C1(2, 1) = (k13)*(k21 - L(1, 0));
+    C1(2, 2) = (k13)*(k21 - L(2, 0));
+
+    C2(2, 0) = C2(2, 1) = C2(2, 2) = (k21)*(k13);
+
+    C3(2, 0) = (k10 + k12 + k13 - L(0, 0))*(k21 - L(0, 0)) - (k21)*(k12);
+    C3(2, 1) = (k10 + k12 + k13 - L(1, 0))*(k21 - L(1, 0)) - (k21)*(k12);
+    C3(2, 2) = (k10 + k12 + k13 - L(2, 0))*(k21 - L(2, 0)) - (k21)*(k12);
+
+    C1(0, 0) = C1(0, 0)/D(0, 0);
+    C1(1, 0) = C1(1, 0)/D(0, 0);
+    C1(2, 0) = C1(2, 0)/D(0, 0);
+
+    C2(0, 0) = C2(0, 0)/D(0, 0);
+    C2(1, 0) = C2(1, 0)/D(0, 0);
+    C2(2, 0) = C2(2, 0)/D(0, 0);
+
+    C3(0, 0) = C3(0, 0)/D(0, 0);
+    C3(1, 0) = C3(1, 0)/D(0, 0);
+    C3(2, 0) = C3(2, 0)/D(0, 0);
+
+    C1(0, 1) = C1(0, 1)/D(1, 0);
+    C1(1, 1) = C1(1, 1)/D(1, 0);
+    C1(2, 1) = C1(2, 1)/D(1, 0);
+
+    C2(0, 1) = C2(0, 1)/D(1, 0);
+    C2(1, 1) = C2(1, 1)/D(1, 0);
+    C2(2, 1) = C2(2, 1)/D(1, 0);
+
+    C3(0, 1) = C3(0, 1)/D(1, 0);
+    C3(1, 1) = C3(1, 1)/D(1, 0);
+    C3(2, 1) = C3(2, 1)/D(1, 0);
+
+    C1(0, 2) = C1(0, 2)/D(2, 0);
+    C1(1, 2) = C1(1, 2)/D(2, 0);
+    C1(2, 2) = C1(2, 2)/D(2, 0);
+
+    C2(0, 2) = C2(0, 2)/D(2, 0);
+    C2(1, 2) = C2(1, 2)/D(2, 0);
+    C2(2, 2) = C2(2, 2)/D(2, 0);
+
+    C3(0, 2) = C3(0, 2)/D(2, 0);
+    C3(1, 2) = C3(1, 2)/D(2, 0);
+    C3(2, 2) = C3(2, 2)/D(2, 0);
+
+    return 1;
+  }
+
+  // now write R interfaces to check the solutions
+}
+
+using namespace Rcpp;
+extern "C" SEXP _rxode2_solComp2cpp(SEXP k10s, SEXP k12s, SEXP k21s) {
+  double k10d = as<double>(k10s);
+  double k12d = as<double>(k12s);
+  double k21d = as<double>(k21s);
+  Eigen::Matrix<double, Eigen::Dynamic, 2> g(2, 2);
+  v = 1.0;
+  k10 = k10d;
+  k12 = k12d;
+  k21 = k21d;
+  Eigen::Matrix<double, 2, 1> L;
+  Eigen::Matrix<double, 2, 2> C1;
+  Eigen::Matrix<double, 2, 2> C2;
+  if (stan::solComp2Cpp(g, L, C1, C2) == 1) {
+    return Rcpp::List::create(Rcpp::_["L"]=Rcpp::wrap(L),
+                              Rcpp::_["C1"]=Rcpp::wrap(C1),
+                              Rcpp::_["C2"]=Rcpp::wrap(C2));
+  }
+  return R_NilValue;
+}
+
+using namespace Rcpp;
+extern "C" SEXP _rxode2_solComp3cpp(SEXP sK10, SEXP sK12, SEXP sK21,
+                                    SEXP sK13, SEXP sK31) {
+  double k10d = as<double>(sK10);
+  double k12d = as<double>(sK12);
+  double k21d = as<double>(sK21);
+  double k13d = as<double>(sK13);
+  double k31d = as<double>(sK31);
+  Eigen::Matrix<double, Eigen::Dynamic, 2> g(3, 2);
+  v = 1.0;
+  k10 = k10d;
+  k12 = k12d;
+  k21 = k21d;
+  k13 = k13d;
+  k31 = k31d;
+  Eigen::Matrix<double, 3, 1> L;
+  Eigen::Matrix<double, 3, 3> C1;
+  Eigen::Matrix<double, 3, 3> C2;
+  Eigen::Matrix<double, 3, 3> C3;
+  if (stan::solComp3Cpp(g, L, C1, C2, C3) == 1) {
+    return Rcpp::List::create(Rcpp::_["L"]=Rcpp::wrap(L),
+                              Rcpp::_["C1"]=Rcpp::wrap(C1),
+                              Rcpp::_["C2"]=Rcpp::wrap(C2),
+                              Rcpp::_["C3"]=Rcpp::wrap(C3));
+  }
+  return R_NilValue;
+}
diff --git a/src/init.c b/src/init.c
index 1d1e3e348..a60db11fa 100644
--- a/src/init.c
+++ b/src/init.c
@@ -364,8 +364,22 @@ SEXP _rxode2_getClassicEvid(SEXP, SEXP, SEXP, SEXP, SEXP,
                             SEXP, SEXP);
 SEXP _rxode2_rxQs(SEXP);
 
+SEXP _rxode2_solComp2(SEXP, SEXP, SEXP);
+SEXP _rxode2_solComp3(SEXP, SEXP, SEXP, SEXP, SEXP);
+
+SEXP _rxode2_comp1c(SEXP, SEXP, SEXP, SEXP, SEXP,
+                    SEXP);
+
+SEXP _rxode2_solComp2cpp(SEXP, SEXP, SEXP);
+SEXP _rxode2_solComp3cpp(SEXP, SEXP, SEXP, SEXP, SEXP);
+
 void R_init_rxode2(DllInfo *info){
   R_CallMethodDef callMethods[]  = {
+    {"_rxode2_solComp2cpp", (DL_FUNC) &_rxode2_solComp2cpp, 3},
+    {"_rxode2_solComp3cpp", (DL_FUNC) &_rxode2_solComp3cpp, 5},
+    {"_rxode2_comp1c", (DL_FUNC) &_rxode2_comp1c, 6},
+    {"_rxode2_solComp2", (DL_FUNC) &_rxode2_solComp2, 3},
+    {"_rxode2_solComp3", (DL_FUNC) &_rxode2_solComp3, 5},
     {"_rxode2_rxode2parseSetRstudio", (DL_FUNC) &_rxode2_rxode2parseSetRstudio, 1},
     {"_rxode2_rxQs", (DL_FUNC) &_rxode2_rxQs, 1},
     {"_rxode2_rxQr", (DL_FUNC) &_rxode2_rxQr, 1},
diff --git a/src/lincmt.c b/src/lincmt.c
index fda4537d3..8fd61ba26 100644
--- a/src/lincmt.c
+++ b/src/lincmt.c
@@ -7,11 +7,9 @@
 #include <Rmath.h>
 #include <R_ext/Rdynload.h>
 
-#define max2( a , b )  ( (a) > (b) ? (a) : (b) )
-#define isSameTime(xout, xp) (fabs((xout)-(xp))  <= DBL_EPSILON*max2(fabs(xout),fabs(xp)))
-
 void _rxode2parse_unprotect(void);
 
+#include "../inst/include/rxode2.h"
 #include "../inst/include/rxode2parse.h"
 
 extern rx_solve rx_global;
@@ -1362,167 +1360,7 @@ static inline void doAdvan(double *A,// Amounts
   }
 }
 
-static inline int parTrans(int *trans,
-                           double *p1, double *v1,
-                           double *p2, double *p3,
-                           double *p4, double *p5,
-                           unsigned int *ncmt, double *rx_k, double *rx_v, double *rx_k12,
-                           double *rx_k21, double *rx_k13, double *rx_k31){
-  double btemp, ctemp, dtemp;
-  if ((*p5) > 0.) {
-    (*ncmt) = 3;
-    switch (*trans) {
-    case 1: // cl v q vp
-      (*rx_k) = (*p1)/(*v1); // k = CL/V
-      (*rx_v) = (*v1);
-      (*rx_k12) = (*p2)/(*v1); // k12 = Q/V
-      (*rx_k21) = (*p2)/(*p3); // k21 = Q/Vp
-      (*rx_k13) = (*p4)/(*v1); // k31 = Q2/V
-      (*rx_k31) = (*p4)/(*p5); // k31 = Q2/Vp2
-      break;
-    case 2: // k=(*p1) v=(*v1) k12=(*p2) k21=(*p3) k13=(*p4) k31=(*p5)
-      (*rx_k) = (*p1);
-      (*rx_v) = (*v1);
-      (*rx_k12) = (*p2);
-      (*rx_k21) = (*p3);
-      (*rx_k13) = (*p4);
-      (*rx_k31) = (*p5);
-      break;
-    case 11:
-#undef beta
-#define A (1/(*v1))
-#define B (*p3)
-#define C (*p5)
-#define alpha (*p1)
-#define beta (*p2)
-#define gamma (*p4)
-      (*ncmt)=3;
-      (*rx_v)=1/(A+B+C);
-      btemp = -(alpha*C + alpha*B + gamma*A + gamma*B + beta*A + beta*C)*(*rx_v);
-      ctemp = (alpha*beta*C + alpha*gamma*B + beta*gamma*A)*(*rx_v);
-      dtemp = sqrt(btemp*btemp-4*ctemp);
-      (*rx_k21) = 0.5*(-btemp+dtemp);
-      (*rx_k31) = 0.5*(-btemp-dtemp);
-      (*rx_k)   = alpha*beta*gamma/(*rx_k21)/(*rx_k31);
-      (*rx_k12) = ((beta*gamma + alpha*beta + alpha*gamma) -
-                   (*rx_k21)*(alpha+beta+gamma) - (*rx_k) * (*rx_k31) + (*rx_k21)*(*rx_k21))/((*rx_k31) - (*rx_k21));
-      (*rx_k13) = alpha + beta + gamma - ((*rx_k) + (*rx_k12) + (*rx_k21) + (*rx_k31));
-      break;
-    case 10:
-#undef A
-#define A (*v1)
-      (*ncmt)=3;
-      (*rx_v)=1/(A+B+C);
-      btemp = -(alpha*C + alpha*B + gamma*A + gamma*B + beta*A + beta*C)*(*rx_v);
-      ctemp = (alpha*beta*C + alpha*gamma*B + beta*gamma*A)*(*rx_v);
-      dtemp = sqrt(btemp*btemp-4*ctemp);
-      (*rx_k21) = 0.5*(-btemp+dtemp);
-      (*rx_k31) = 0.5*(-btemp-dtemp);
-      (*rx_k)   = alpha*beta*gamma/(*rx_k21)/(*rx_k31);
-      (*rx_k12) = ((beta*gamma + alpha*beta + alpha*gamma) -
-                   (*rx_k21)*(alpha+beta+gamma) - (*rx_k) * (*rx_k31) + (*rx_k21)*(*rx_k21))/((*rx_k31) - (*rx_k21));
-      (*rx_k13) = alpha + beta + gamma - ((*rx_k) + (*rx_k12) + (*rx_k21) + (*rx_k31));
-#undef A
-#undef B
-#undef C
-#undef alpha
-#undef beta
-#undef gamma
-#define beta Rf_beta
-      break;
-    default:
-      return NA_REAL;
-    }
-  } else if ((*p3) > 0.) {
-    (*ncmt) = 2;
-    switch (*trans){
-    case 1: // cl=(*p1) v=(*v1) q=(*p2) vp=(*p3)
-      (*rx_k) = (*p1)/(*v1); // k = CL/V
-      (*rx_v) = (*v1);
-      (*rx_k12) = (*p2)/(*v1); // k12 = Q/V
-      (*rx_k21) = (*p2)/(*p3); // k21 = Q/Vp
-      break;
-    case 2: // k=(*p1), (*v1)=v k12=(*p2) k21=(*p3)
-      (*rx_k) = (*p1);
-      (*rx_v) = (*v1);
-      (*rx_k12) = (*p2);
-      (*rx_k21) = (*p3);
-      break;
-    case 3: // cl=(*p1) v=(*v1) q=(*p2) vss=(*p3)
-      (*rx_k) = (*p1)/(*v1); // k = CL/V
-      (*rx_v) = (*v1);
-      (*rx_k12) = (*p2)/(*v1); // k12 = Q/V
-      (*rx_k21) = (*p2)/((*p3)-(*v1)); // k21 = Q/(Vss-V)
-      break;
-    case 4: // alpha=(*p1) beta=(*p2) k21=(*p3)
-      (*rx_v) = (*v1);
-      (*rx_k21) = (*p3);
-      (*rx_k) = (*p1)*(*p2)/(*rx_k21); // (*p1) = alpha (*p2) = beta
-      (*rx_k12) = (*p1) + (*p2) - (*rx_k21) - (*rx_k);
-      break;
-    case 5: // alpha=(*p1) beta=(*p2) aob=(*p3)
-      (*rx_v)=(*v1);
-      (*rx_k21) = ((*p3)*(*p2)+(*p1))/((*p3)+1.0);
-      (*rx_k) = ((*p1)*(*p2))/(*rx_k21);
-      (*rx_k12) = (*p1) + (*p2) - (*rx_k21) - (*rx_k);
-      break;
-    case 11: // A2 V, alpha=(*p1), beta=(*p2), k21
-#undef beta
-#define A (1/(*v1))
-#define B (*p3)
-#define alpha (*p1)
-#define beta (*p2)
-      (*ncmt)=2;
-      (*rx_v)   = 1/(A+B);
-      (*rx_k21) = (A*beta + B*alpha)*(*rx_v);
-      (*rx_k)   = alpha*beta/(*rx_k21);
-      (*rx_k12) = alpha+beta-(*rx_k21)-(*rx_k);
-      break;
-    case 10: // A=(*v1), alpha=(*p1), beta=(*p2), B=(*p3)
-      // Convert to A (right now A=(*v1) or A=1/(*v1))
-#undef A
-#define A (*v1)
-      (*ncmt)=2;
-      (*rx_v)   = 1/(A + B);
-      (*rx_k21) = (A*beta + B*alpha)*(*rx_v);
-      (*rx_k)   = alpha*beta/(*rx_k21);
-      (*rx_k12) = alpha + beta - (*rx_k21) - (*rx_k);
-#undef A
-#undef B
-#undef alpha
-#undef beta
-#define beta Rf_beta
-      break;
-    default:
-      return NA_REAL;
-    }
-  } else if ((*p1) > 0.) {
-    (*ncmt) = 1;
-    switch(*trans){
-    case 1: // cl v
-      (*rx_k) = (*p1)/(*v1); // k = CL/V
-      (*rx_v) = (*v1);
-      break;
-    case 2: // k V
-      (*rx_k) = (*p1);
-      (*rx_v) = (*v1);
-      break;
-    case 11: // alpha V
-      (*rx_k) = (*p1);
-      (*rx_v) = (*v1);
-      break;
-    case 10: // alpha A
-      (*rx_k) = (*p1);
-      (*rx_v) = 1/(*v1);
-      break;
-    default:
-      return 0;
-    }
-  } else {
-    return 0;
-  }
-  return 1;
-}
+#include "parTrans.h"
 
 void linCmtPar1(double *v, double *k,
                 double *vss,
diff --git a/src/parTrans.h b/src/parTrans.h
new file mode 100644
index 000000000..f746fa1b2
--- /dev/null
+++ b/src/parTrans.h
@@ -0,0 +1,161 @@
+static inline int parTrans(int *trans,
+                           double *p1, double *v1,
+                           double *p2, double *p3,
+                           double *p4, double *p5,
+                           unsigned int *ncmt, double *rx_k, double *rx_v, double *rx_k12,
+                           double *rx_k21, double *rx_k13, double *rx_k31){
+  double btemp, ctemp, dtemp;
+  if ((*p5) > 0.) {
+    (*ncmt) = 3;
+    switch (*trans) {
+    case 1: // cl v q vp
+      (*rx_k) = (*p1)/(*v1); // k = CL/V
+      (*rx_v) = (*v1);
+      (*rx_k12) = (*p2)/(*v1); // k12 = Q/V
+      (*rx_k21) = (*p2)/(*p3); // k21 = Q/Vp
+      (*rx_k13) = (*p4)/(*v1); // k31 = Q2/V
+      (*rx_k31) = (*p4)/(*p5); // k31 = Q2/Vp2
+      break;
+    case 2: // k=(*p1) v=(*v1) k12=(*p2) k21=(*p3) k13=(*p4) k31=(*p5)
+      (*rx_k) = (*p1);
+      (*rx_v) = (*v1);
+      (*rx_k12) = (*p2);
+      (*rx_k21) = (*p3);
+      (*rx_k13) = (*p4);
+      (*rx_k31) = (*p5);
+      break;
+    case 11:
+#undef beta
+#define A (1/(*v1))
+#define B (*p3)
+#define C (*p5)
+#define alpha (*p1)
+#define beta (*p2)
+#define gamma (*p4)
+      (*ncmt)=3;
+      (*rx_v)=1/(A+B+C);
+      btemp = -(alpha*C + alpha*B + gamma*A + gamma*B + beta*A + beta*C)*(*rx_v);
+      ctemp = (alpha*beta*C + alpha*gamma*B + beta*gamma*A)*(*rx_v);
+      dtemp = sqrt(btemp*btemp-4*ctemp);
+      (*rx_k21) = 0.5*(-btemp+dtemp);
+      (*rx_k31) = 0.5*(-btemp-dtemp);
+      (*rx_k)   = alpha*beta*gamma/(*rx_k21)/(*rx_k31);
+      (*rx_k12) = ((beta*gamma + alpha*beta + alpha*gamma) -
+                   (*rx_k21)*(alpha+beta+gamma) - (*rx_k) * (*rx_k31) + (*rx_k21)*(*rx_k21))/((*rx_k31) - (*rx_k21));
+      (*rx_k13) = alpha + beta + gamma - ((*rx_k) + (*rx_k12) + (*rx_k21) + (*rx_k31));
+      break;
+    case 10:
+#undef A
+#define A (*v1)
+      (*ncmt)=3;
+      (*rx_v)=1/(A+B+C);
+      btemp = -(alpha*C + alpha*B + gamma*A + gamma*B + beta*A + beta*C)*(*rx_v);
+      ctemp = (alpha*beta*C + alpha*gamma*B + beta*gamma*A)*(*rx_v);
+      dtemp = sqrt(btemp*btemp-4*ctemp);
+      (*rx_k21) = 0.5*(-btemp+dtemp);
+      (*rx_k31) = 0.5*(-btemp-dtemp);
+      (*rx_k)   = alpha*beta*gamma/(*rx_k21)/(*rx_k31);
+      (*rx_k12) = ((beta*gamma + alpha*beta + alpha*gamma) -
+                   (*rx_k21)*(alpha+beta+gamma) - (*rx_k) * (*rx_k31) + (*rx_k21)*(*rx_k21))/((*rx_k31) - (*rx_k21));
+      (*rx_k13) = alpha + beta + gamma - ((*rx_k) + (*rx_k12) + (*rx_k21) + (*rx_k31));
+#undef A
+#undef B
+#undef C
+#undef alpha
+#undef beta
+#undef gamma
+#define beta Rf_beta
+      break;
+    default:
+      return NA_REAL;
+    }
+  } else if ((*p3) > 0.) {
+    (*ncmt) = 2;
+    switch (*trans){
+    case 1: // cl=(*p1) v=(*v1) q=(*p2) vp=(*p3)
+      (*rx_k) = (*p1)/(*v1); // k = CL/V
+      (*rx_v) = (*v1);
+      (*rx_k12) = (*p2)/(*v1); // k12 = Q/V
+      (*rx_k21) = (*p2)/(*p3); // k21 = Q/Vp
+      break;
+    case 2: // k=(*p1), (*v1)=v k12=(*p2) k21=(*p3)
+      (*rx_k) = (*p1);
+      (*rx_v) = (*v1);
+      (*rx_k12) = (*p2);
+      (*rx_k21) = (*p3);
+      break;
+    case 3: // cl=(*p1) v=(*v1) q=(*p2) vss=(*p3)
+      (*rx_k) = (*p1)/(*v1); // k = CL/V
+      (*rx_v) = (*v1);
+      (*rx_k12) = (*p2)/(*v1); // k12 = Q/V
+      (*rx_k21) = (*p2)/((*p3)-(*v1)); // k21 = Q/(Vss-V)
+      break;
+    case 4: // alpha=(*p1) beta=(*p2) k21=(*p3)
+      (*rx_v) = (*v1);
+      (*rx_k21) = (*p3);
+      (*rx_k) = (*p1)*(*p2)/(*rx_k21); // (*p1) = alpha (*p2) = beta
+      (*rx_k12) = (*p1) + (*p2) - (*rx_k21) - (*rx_k);
+      break;
+    case 5: // alpha=(*p1) beta=(*p2) aob=(*p3)
+      (*rx_v)=(*v1);
+      (*rx_k21) = ((*p3)*(*p2)+(*p1))/((*p3)+1.0);
+      (*rx_k) = ((*p1)*(*p2))/(*rx_k21);
+      (*rx_k12) = (*p1) + (*p2) - (*rx_k21) - (*rx_k);
+      break;
+    case 11: // A2 V, alpha=(*p1), beta=(*p2), k21
+#undef beta
+#define A (1/(*v1))
+#define B (*p3)
+#define alpha (*p1)
+#define beta (*p2)
+      (*ncmt)=2;
+      (*rx_v)   = 1/(A+B);
+      (*rx_k21) = (A*beta + B*alpha)*(*rx_v);
+      (*rx_k)   = alpha*beta/(*rx_k21);
+      (*rx_k12) = alpha+beta-(*rx_k21)-(*rx_k);
+      break;
+    case 10: // A=(*v1), alpha=(*p1), beta=(*p2), B=(*p3)
+      // Convert to A (right now A=(*v1) or A=1/(*v1))
+#undef A
+#define A (*v1)
+      (*ncmt)=2;
+      (*rx_v)   = 1/(A + B);
+      (*rx_k21) = (A*beta + B*alpha)*(*rx_v);
+      (*rx_k)   = alpha*beta/(*rx_k21);
+      (*rx_k12) = alpha + beta - (*rx_k21) - (*rx_k);
+#undef A
+#undef B
+#undef alpha
+#undef beta
+#define beta Rf_beta
+      break;
+    default:
+      return NA_REAL;
+    }
+  } else if ((*p1) > 0.) {
+    (*ncmt) = 1;
+    switch(*trans){
+    case 1: // cl v
+      (*rx_k) = (*p1)/(*v1); // k = CL/V
+      (*rx_v) = (*v1);
+      break;
+    case 2: // k V
+      (*rx_k) = (*p1);
+      (*rx_v) = (*v1);
+      break;
+    case 11: // alpha V
+      (*rx_k) = (*p1);
+      (*rx_v) = (*v1);
+      break;
+    case 10: // alpha A
+      (*rx_k) = (*p1);
+      (*rx_v) = 1/(*v1);
+      break;
+    default:
+      return 0;
+    }
+  } else {
+    return 0;
+  }
+  return 1;
+}
diff --git a/src/par_solve.cpp b/src/par_solve.cpp
index c530fc7b6..3cabb1914 100644
--- a/src/par_solve.cpp
+++ b/src/par_solve.cpp
@@ -14,8 +14,6 @@
 #include "../inst/include/rxode2parseHandleEvid.h"
 #include "../inst/include/rxode2parseGetTime.h"
 #define SORT gfx::timsort
-#define isSameTimeOp(xout, xp) (op->stiff == 0 ? isSameTimeDop(xout, xp) : isSameTime(xout, xp))
-// dop853 is same time
 
 extern "C" {
 #include "dop853.h"
diff --git a/src/solComp.c b/src/solComp.c
new file mode 100644
index 000000000..db4f7aad7
--- /dev/null
+++ b/src/solComp.c
@@ -0,0 +1,80 @@
+#define USE_FC_LEN_T
+#define STRICT_R_HEADERS
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>   /* dj: import intptr_t */
+//#include "ode.h"
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
+#include <Rmath.h>
+#ifdef ENABLE_NLS
+#include <libintl.h>
+#define _(String) dgettext ("rxode2", String)
+/* replace pkg as appropriate */
+#else
+#define _(String) (String)
+#endif
+#include "solComp.h"
+
+SEXP _rxode2_solComp2(SEXP sK10, SEXP sK12, SEXP sK21) {
+  int pro = 0;
+  SEXP L = PROTECT(allocVector(REALSXP, 2)); pro++;
+  SEXP C1 = PROTECT(allocVector(REALSXP, 4)); pro++;
+  SEXP dm = PROTECT(allocVector(INTSXP, 2)); pro++;
+  SEXP C2 = PROTECT(allocVector(REALSXP, 4)); pro++;
+  int *dmi = INTEGER(dm);
+  dmi[0] = dmi[1] = 2;
+  Rf_setAttrib(C1, R_DimSymbol, dm);
+  Rf_setAttrib(C2, R_DimSymbol, dm);
+  if (!solComp2C(REAL(sK10), REAL(sK12), REAL(sK21),
+                REAL(L), REAL(C1), REAL(C2))) {
+    UNPROTECT(pro);
+    return R_NilValue;
+  }
+  SEXP lst   = PROTECT(allocVector(VECSXP, 3)); pro++;
+  SEXP names = PROTECT(allocVector(STRSXP, 3)); pro++;
+  SET_STRING_ELT(names,0,mkChar("L"));
+  SET_VECTOR_ELT(lst,  0, L);
+  SET_STRING_ELT(names,1,mkChar("C1"));
+  SET_VECTOR_ELT(lst,  1, C1);
+  SET_STRING_ELT(names,2,mkChar("C2"));
+  SET_VECTOR_ELT(lst,  2, C2);
+  Rf_setAttrib(lst, R_NamesSymbol, names);
+  UNPROTECT(pro);
+  return lst;
+}
+
+SEXP _rxode2_solComp3(SEXP sK10, SEXP sK12, SEXP sK21, SEXP sK13, SEXP sK31) {
+  int pro=0;
+  SEXP L = PROTECT(allocVector(REALSXP, 3)); pro++;
+  SEXP C1 = PROTECT(allocVector(REALSXP, 9)); pro++;
+  SEXP dm = PROTECT(allocVector(INTSXP, 2)); pro++;
+  SEXP C2 = PROTECT(allocVector(REALSXP, 9)); pro++;
+  SEXP C3 = PROTECT(allocVector(REALSXP, 9)); pro++;
+  int *dmi = INTEGER(dm);
+  dmi[0] = dmi[1] = 3;
+  Rf_setAttrib(C1, R_DimSymbol, dm);
+  Rf_setAttrib(C2, R_DimSymbol, dm);
+  Rf_setAttrib(C3, R_DimSymbol, dm);
+  if (!solComp3C(REAL(sK10), REAL(sK12), REAL(sK21), REAL(sK13), REAL(sK31),
+                 REAL(L), REAL(C1), REAL(C2), REAL(C3))) {
+    UNPROTECT(pro);
+    return R_NilValue;
+  }
+  SEXP lst   = PROTECT(allocVector(VECSXP, 4)); pro++;
+  SEXP names = PROTECT(allocVector(STRSXP, 4)); pro++;
+  SET_STRING_ELT(names,0,mkChar("L"));
+  SET_VECTOR_ELT(lst,  0, L);
+  SET_STRING_ELT(names,1,mkChar("C1"));
+  SET_VECTOR_ELT(lst,  1, C1);
+  SET_STRING_ELT(names,2,mkChar("C2"));
+  SET_VECTOR_ELT(lst,  2, C2);
+  SET_STRING_ELT(names,3,mkChar("C3"));
+  SET_VECTOR_ELT(lst,  3, C3);
+  Rf_setAttrib(lst, R_NamesSymbol, names);
+  UNPROTECT(pro);
+  return lst;
+}
diff --git a/src/solComp.h b/src/solComp.h
new file mode 100644
index 000000000..f9d687ce7
--- /dev/null
+++ b/src/solComp.h
@@ -0,0 +1,128 @@
+#ifndef __SOLCOMP_H__
+#define __SOLCOMP_H__
+
+static inline int solComp2C(double *k10, double *k12, double *k21,
+                            double *L, double *C1, double *C2) {
+  // blas and R uses row major matrices
+  double sum = (*k10) + (*k12) + (*k21);
+  double disc= sqrt(sum*sum - 4.0* (*k10)*(*k21));
+  double div[2];
+  double tmp;
+  L[0] = 0.5*(sum + disc);
+  L[1] = 0.5*(sum - disc);
+  div[0] = L[1] - L[0];
+  div[1] = L[0] - L[1];
+  if (div[0]*div[1] == 0) return 0;
+  C1[0] = (*k21) - L[0];
+  C1[2] = (*k21) - L[1];
+  C2[0] = C2[2] = (*k21);
+  C1[1] = C1[3] = (*k12);
+  tmp = (*k10) + (*k12);
+  C2[1] = tmp - L[0];
+  C2[3] = tmp - L[1];
+  C1[0] = C1[0]/div[0];
+  C1[1] = C1[1]/div[0];
+  C2[0] = C2[0]/div[0];
+  C2[1] = C2[1]/div[0];
+  C1[2] = C1[2]/div[1];
+  C1[3] = C1[3]/div[1];
+  C2[2] = C2[2]/div[1];
+  C2[3] = C2[3]/div[1];
+  return 1;
+}
+
+static inline int solComp3C(double *k10, double *k12, double *k21,
+                            double *k13, double *k31,
+                            double *L, double *C1, double *C2, double *C3) {
+  // Get Lambdas
+  double A1 = *k10 + *k12 + *k13 + *k21 + *k31;
+  double A2 = (*k10)* (*k21) + (*k10)*(*k31) +
+    (*k12)*(*k31) + (*k13)*(*k21) + (*k21)*(*k31);
+  double A3 = (*k21)*(*k31)*(*k10);
+  double Q  = (A1*A1 - 3.0*A2)/9.0;
+  double RQ = 2.0*sqrt(Q);
+  double R  = (2.0*A1*A1*A1 - 9.0*A1*A2 + 27.0*A3)/54.0;
+  double M  = Q*Q*Q - R*R;
+  if (M < 0) return 0;//stop("Error: Not real roots.")
+  double Th = acos(8.0*R/(RQ*RQ*RQ));
+  L[0] = RQ*cos(Th/3.0) + A1/3.0;
+  L[1] = RQ*cos((Th + M_2PI)/3.0) + A1/3.0;
+  L[2] = RQ*cos((Th + 2*M_2PI)/3.0) + A1/3.0;
+  double D[3];
+  D[0] = (L[1] - L[0])*(L[2] - L[0]);
+  D[1] = (L[0] - L[1])*(L[2] - L[1]);
+  D[2] = (L[0] - L[2])*(L[1] - L[2]);
+  if (D[0]*D[1]*D[2] == 0.0) return 0;
+
+  C1[0] = (*k21 - L[0])*(*k31 - L[0]);
+  C1[3] = (*k21 - L[1])*(*k31 - L[1]);
+  C1[6] = (*k21 - L[2])*(*k31 - L[2]);
+
+  C2[0] = (*k21)*(*k31 - L[0]);
+  C2[3] = (*k21)*(*k31 - L[1]);
+  C2[6] = (*k21)*(*k31 - L[2]);
+
+  C3[0] = (*k31)*(*k21 - L[0]);
+  C3[3] = (*k31)*(*k21 - L[1]);
+  C3[6] = (*k31)*(*k21 - L[2]);
+
+  C1[1] = (*k12)*(*k31 - L[0]);
+  C1[4] = (*k12)*(*k31 - L[1]);
+  C1[7] = (*k12)*(*k31 - L[2]);
+
+  C2[1] = (*k10 + *k12 + *k13 - L[0])*(*k31 - L[0]) - (*k31)*(*k13);
+  C2[4] = (*k10 + *k12 + *k13 - L[1])*(*k31 - L[1]) - (*k31)*(*k13);
+  C2[7] = (*k10 + *k12 + *k13 - L[2])*(*k31 - L[2]) - (*k31)*(*k13);
+
+  C3[1] = C3[4] = C3[7] = (*k12)*(*k31);
+
+  C1[2] = (*k13)*(*k21 - L[0]);
+  C1[5] = (*k13)*(*k21 - L[1]);
+  C1[8] = (*k13)*(*k21 - L[2]);
+
+  C2[2] = C2[5] = C2[8] = (*k21)*(*k13);
+
+  C3[2] = (*k10 + *k12 + *k13 - L[0])*(*k21 - L[0]) - (*k21)*(*k12);
+  C3[5] = (*k10 + *k12 + *k13 - L[1])*(*k21 - L[1]) - (*k21)*(*k12);
+  C3[8] = (*k10 + *k12 + *k13 - L[2])*(*k21 - L[2]) - (*k21)*(*k12);
+
+  C1[0] = C1[0]/D[0];
+  C1[1] = C1[1]/D[0];
+  C1[2] = C1[2]/D[0];
+
+  C2[0] = C2[0]/D[0];
+  C2[1] = C2[1]/D[0];
+  C2[2] = C2[2]/D[0];
+
+  C3[0] = C3[0]/D[0];
+  C3[1] = C3[1]/D[0];
+  C3[2] = C3[2]/D[0];
+
+  C1[3] = C1[3]/D[1];
+  C1[4] = C1[4]/D[1];
+  C1[5] = C1[5]/D[1];
+
+  C2[3] = C2[3]/D[1];
+  C2[4] = C2[4]/D[1];
+  C2[5] = C2[5]/D[1];
+
+  C3[3] = C3[3]/D[1];
+  C3[4] = C3[4]/D[1];
+  C3[5] = C3[5]/D[1];
+
+  C1[6] = C1[6]/D[2];
+  C1[7] = C1[7]/D[2];
+  C1[8] = C1[8]/D[2];
+
+  C2[6] = C2[6]/D[2];
+  C2[7] = C2[7]/D[2];
+  C2[8] = C2[8]/D[2];
+
+  C3[6] = C3[6]/D[2];
+  C3[7] = C3[7]/D[2];
+  C3[8] = C3[8]/D[2];
+
+  return 1;
+}
+
+#endif
diff --git a/tests/testthat/test-solComp.R b/tests/testthat/test-solComp.R
new file mode 100644
index 000000000..d0f82fe26
--- /dev/null
+++ b/tests/testthat/test-solComp.R
@@ -0,0 +1,35 @@
+test_that("test the matrices of the linear compartment solutions .solComp2", {
+
+  v <- list(L = c(4.07546291005291, 0.0245370899470903),
+            C1 = structure(c(0.759200006770939,-0.740571447916969,
+                             0.240799993229061, 0.740571447916969),
+                           dim = c(2L, 2L)),
+            C2 = structure(c(-0.246857149305656, 0.240799993229061,
+                             0.246857149305656, 0.759200006770939),
+                           dim = c(2L, 2L)))
+
+  expect_equal(.solComp2(k10=0.1, k12=3, k21=1, cpp=FALSE), v)
+
+  expect_equal(.solComp2(k10=0.1, k12=3, k21=1, cpp=TRUE), v)
+
+  .Call(`_rxode2_solComp2cpp`, 0.1, 3, 1)
+
+  v <- list(L = c(5.8977832399092, 0.0122878774411843, 0.689928882649618),
+            C1 = structure(c(0.86252789433926, -0.528317313419071, -0.319585969277925,
+                             0.120784508040635, 0.366861472939276, 0.495310665672737,
+                             0.016687597620105, 0.161455840479795, -0.175724696394812),
+                           dim = c(3L, 3L)),
+            C2 = structure(c(-0.17610577113969, 0.107868659665074, 0.0652511460029109,
+                             0.122287157646425, 0.371425504011093, 0.501472700759773,
+                             0.053818613493265, 0.520705836323833, -0.566723846762684),
+                           dim = c(3L, 3L)),
+            C3 = structure(c(-0.0798964923194813, 0.0489383595021832, 0.0296034459956658,
+                             0.123827666418184, 0.37610452556983, 0.507789987948271,
+                             -0.043931174098703, -0.425042885072013, 0.462606566056063),
+                           dim = c(3L,  3L)))
+
+  expect_equal(.solComp3(k10=0.1, k12=3, k21=1, k13=2, k31=0.5, cpp=FALSE), v)
+
+  expect_equal(.solComp3(k10=0.1, k12=3, k21=1, k13=2, k31=0.5, cpp=TRUE), v)
+
+})