diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 1271156..606d6eb 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -239,8 +239,8 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
     A_matrix.insert(1, 1) = 1;
   } else {
     double sx = (alpha[0] * time_step) / (dx * dx);
-    A_matrix.insert(1, 1) = -1. - 2. * sx;
-    A_matrix.insert(1, 0) = sx;
+    A_matrix.insert(1, 1) = -1. - 3. * sx;
+    A_matrix.insert(1, 0) = 2. * sx;
     A_matrix.insert(1, 2) = sx;
   }
 
@@ -261,9 +261,9 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
     A_matrix.insert(A_size - 2, A_size - 2) = 1;
   } else {
     double sx = (alpha[size - 1] * time_step) / (dx * dx);
-    A_matrix.insert(A_size - 2, A_size - 2) = -1. - 2. * sx;
+    A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
     A_matrix.insert(A_size - 2, A_size - 3) = sx;
-    A_matrix.insert(A_size - 2, A_size - 1) = sx;
+    A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
   }
 
   A_matrix.insert(A_size - 1, A_size - 1) = 1;
@@ -294,15 +294,14 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow(
     b_vector[j + 1] = -c[j] - t0_c_j;
   }
 
-  if (!left_constant) {
-    // this is not correct currently.We will fix this when we are able to define
-    // FLUX boundary conditions
-    b_vector[0] = getBCFromFlux(left, c[0], alpha[0]);
-  }
+  // this is not correct currently.We will fix this when we are able to define
+  // FLUX boundary conditions
+  b_vector[0] =
+      (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
 
-  if (!right_constant) {
-    b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]);
-  }
+  b_vector[b_size - 1] =
+      (right_constant ? right.value
+                      : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
 }
 
 void Diffusion::BTCSDiffusion::setTimestep(double time_step) {