diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index fd1b5a7..af5907d 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -257,27 +257,26 @@ void BTCSDiffusion::fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
         getBCFromFlux(right, c(row, ncol - 1), alpha[ncol - 1]);
   }
 
-  for (int j = 1; j < offset - 1; j++) {
-    boundary_condition tmp_bc = this->bc(row, j-1);
+  for (int j = 0; j < ncol; j++) {
+    boundary_condition tmp_bc = this->bc(row, j);
 
     if (tmp_bc.type == BTCSDiffusion::BC_CONSTANT) {
-      b_vector[offset * row + j] = tmp_bc.value;
+      b_vector[offset * row + (j+1)] = tmp_bc.value;
       continue;
     }
 
     double y_values[3];
     y_values[0] =
-        (row != 0 ? c(row - 1, j - 1)
-                  : getBCFromFlux(tmp_bc, c(row, j - 1), alpha[j - 1]));
-    y_values[1] = c(row, j - 1);
+        (row != 0 ? c(row - 1, j) : getBCFromFlux(tmp_bc, c(row, j), alpha[j]));
+    y_values[1] = c(row, j);
     y_values[2] =
-        (row != nrow - 1 ? c(row + 1, j - 1)
-                         : getBCFromFlux(tmp_bc, c(row, j - 1), alpha[j - 1]));
+        (row != nrow - 1 ? c(row + 1, j)
+                         : getBCFromFlux(tmp_bc, c(row, j), alpha[j]));
 
     double t0_c =
-        alpha[j - 1] *
+        alpha[j] *
         ((y_values[0] - 2 * y_values[1] + y_values[2]) / (delta * delta));
-    b_vector[offset * row + j] = -c(row, j - 1) - t0_c;
+    b_vector[offset * row + (j + 1)] = -c(row, j) - t0_c;
   }
 }