diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 53667d2..f4888ee 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -138,6 +138,50 @@ void BTCSDiffusion::simulate1D(Eigen::Map<Eigen::VectorXd> &c,
   c = x_vector.segment(!left_is_constant, c.size());
 }
 
+void BTCSDiffusion::simulate2D(Eigen::Map<Eigen::MatrixXd> &c,
+                               Eigen::Map<const Eigen::MatrixXd> &alpha) {
+
+  int n_cols = c.cols();
+  unsigned int size = (this->grid_cells[0] + 2) * (this->grid_cells[1]);
+
+  A_matrix.resize(size, size);
+  A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
+
+  for (int i = 0; i < c.rows(); i++) {
+    bool left = bc[i*n_cols].type == BTCSDiffusion::BC_CONSTANT;
+    bool right = bc[((i+1)*n_cols)-1].type == BTCSDiffusion::BC_CONSTANT;
+    fillMatrixFromRow(alpha, i, left, right, domain_size[0]);
+  }
+}
+
+inline void
+BTCSDiffusion::fillMatrixFromRow(Eigen::Map<const Eigen::MatrixXd> &alpha,
+                                 int row, bool left_constant,
+                                 bool right_constant, int delta) {
+
+  int n_cols = A_matrix.cols();
+  int offset = n_cols * row;
+
+  A_matrix.insert(offset, offset) = !left_constant;
+
+  if (left_constant)
+    A_matrix.insert(offset + 1, offset + 1) = 1;
+
+  A_matrix.insert(offset + (n_cols - 1), offset + (n_cols - 1)) =
+      !right_constant;
+
+  if (right_constant)
+    A_matrix.insert(offset + (n_cols - 2), offset + (n_cols - 2)) = 1;
+
+  for (int j = 1 + left_constant; j < offset - (1 - right_constant); j++) {
+    double sx = (alpha(row, j) * this->time_step) / (delta * delta);
+
+    A_matrix.insert(offset + j, offset + j) = -1. - 2. * sx;
+    A_matrix.insert(offset + j, offset + (j - 1)) = sx;
+    A_matrix.insert(offset + j, offset + (j + 1)) = sx;
+  }
+}
+
 void BTCSDiffusion::setTimestep(double time_step) {
   this->time_step = time_step;
 }
@@ -151,7 +195,12 @@ void BTCSDiffusion::simulate(std::vector<double> &c,
                this->grid_cells[0]);
   }
   if (this->grid_dim == 2) {
+    assert(c.size() == grid_cells[0] * grid_cells[1]);
+    Eigen::Map<Eigen::MatrixXd> c_in(c.data(), this->grid_cells[1],
+                                     this->grid_cells[0]);
 
+    Eigen::Map<const Eigen::MatrixXd> alpha_in(
+        alpha.data(), this->grid_cells[1], this->grid_cells[0]);
   }
 }
 
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 072dff8..f6be297 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -140,11 +140,17 @@ private:
   void simulate1D(Eigen::Map<Eigen::VectorXd> &c, boundary_condition left,
                   boundary_condition right, const std::vector<double> &alpha,
                   double dx, int size);
-  void simulate2D(std::vector<double> &c);
+  void simulate2D(Eigen::Map<Eigen::MatrixXd> &c,
+                  Eigen::Map<const Eigen::MatrixXd> &alpha);
+
+  inline void fillMatrixFromRow(Eigen::Map<const Eigen::MatrixXd> &alpha,
+                                int row, bool left_constant,
+                                bool right_constant, int delta);
+
   void simulate3D(std::vector<double> &c);
   inline double getBCFromFlux(boundary_condition bc, double nearest_value,
                               double neighbor_alpha);
-    inline void solveLES();
+  inline void solveLES();
   void updateInternals();
 
   std::vector<boundary_condition> bc;