diff --git a/src/compB.cpp b/src/compB.cpp
index 129318791..96974b4e6 100644
--- a/src/compB.cpp
+++ b/src/compB.cpp
@@ -32,7 +32,7 @@ namespace stan {
     Eigen::Matrix<T, Eigen::Dynamic, 2>
     parTransB(const Eigen::Matrix<T, Eigen::Dynamic, 1>& p,
                   const int& ncmt,
-                  const int& trans){
+                  const int& trans) {
       Eigen::Matrix<T, Eigen::Dynamic, 2> g(ncmt,3);
       T btemp, ctemp, dtemp;
 #define p1    p[0]
@@ -225,11 +225,154 @@ namespace stan {
 #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) {
+
+    T sum = k10 + k12 + k21;
+    T disc= sqrt(sum*sum - 4.0 * k10 * k21);
+
+    Eigen::Matrix<T, 2, 1> div;
+
+    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;
+
+    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;
+
+    T tmp = k10 + k12;
+
+    C2(1, 0) = tmp - L(0, 0);
+    C2(0, 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;
+  }
 }