diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 4df37a9..4f7b109 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -141,11 +141,10 @@ void Diffusion::BTCSDiffusion::simulate2D(
   int n_rows = this->grid_cells[1];
   int n_cols = this->grid_cells[0];
   double dx = this->deltas[0];
-  DMatrixRowMajor t0_c;
 
   double local_dt = this->time_step / BTCS_2D_DT_SIZE;
 
-  t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
+  DMatrixRowMajor t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
 
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_rows; i++) {
@@ -184,12 +183,17 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
 
   // first, iterate over first row
   for (int j = 0; j < n_cols; j++) {
-    y_values[0] = getBCFromFlux(bc(0, j), c(0, j), alpha(0, j));
+    boundary_condition tmp_bc = bc(0,j+1);
+
+    if (tmp_bc.type == Diffusion::BC_CLOSED)
+      continue;
+
+    y_values[0] = getBCFromFlux(tmp_bc, c(0, j), alpha(0, j));
     y_values[1] = c(0, j);
     y_values[2] = c(1, j);
 
     t0_c(0, j) = time_step * alpha(0, j) *
-                 (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
+                 (2*y_values[0] - 3 * y_values[1] + y_values[2]) / (dx * dx);
   }
 
 // then iterate over inlet
@@ -210,12 +214,17 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
 
   // and finally over last row
   for (int j = 0; j < n_cols; j++) {
+    boundary_condition tmp_bc = bc(end+1,j+1);
+
+    if (tmp_bc.type == Diffusion::BC_CLOSED)
+      continue;
+
     y_values[0] = c(end - 1, j);
     y_values[1] = c(end, j);
-    y_values[2] = getBCFromFlux(bc(end, j), c(end, j), alpha(end, j));
+    y_values[2] = getBCFromFlux(tmp_bc, c(end, j), alpha(end, j));
 
     t0_c(end, j) = time_step * alpha(end, j) *
-                   (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
+                   (y_values[0] - 3 * y_values[1] + 2*y_values[2]) / (dx * dx);
   }
 
   return t0_c;
@@ -240,7 +249,7 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
   } else {
     double sx = (alpha[0] * time_step) / (dx * dx);
     A_matrix.insert(1, 1) = -1. - 3. * sx;
-    A_matrix.insert(1, 0) = 2.*sx;
+    A_matrix.insert(1, 0) = 2. * sx;
     A_matrix.insert(1, 2) = sx;
   }
 
@@ -263,7 +272,7 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
     double sx = (alpha[size - 1] * time_step) / (dx * dx);
     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) = 2.*sx;
+    A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
   }
 
   A_matrix.insert(A_size - 1, A_size - 1) = 1;
@@ -294,15 +303,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) {