diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 92855dc..7f121e2 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -129,14 +129,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
     b_vector[i] = -c[j];
   }
 
-  // start to solve
-  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
-      solver;
-  solver.analyzePattern(A_matrix);
-
-  solver.factorize(A_matrix);
-
-  x_vector = solver.solve(b_vector);
+  solveLES();
 
   //fill solution back in place into =c= vector
   for (int i = 0, j = i + !left_is_constant; i < c.size(); i++, j++) {
@@ -154,6 +147,9 @@ void BTCSDiffusion::simulate(std::vector<double> &c,
     simulate1D(c, bc[0], bc[grid_cells[0] + 1], alpha, this->deltas[0],
                this->grid_cells[0]);
   }
+  if (this->grid_dim == 2) {
+
+  }
 }
 
 inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc,
@@ -179,3 +175,14 @@ void BTCSDiffusion::setBoundaryCondition(int index, bctype type, double value) {
   bc[index].type = type;
   bc[index].value = value;
 }
+
+inline void BTCSDiffusion::solveLES() {
+  // start to solve
+  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
+      solver;
+  solver.analyzePattern(A_matrix);
+
+  solver.factorize(A_matrix);
+
+  x_vector = solver.solve(b_vector);
+}
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index ea3473d..88db9f7 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -142,6 +142,7 @@ private:
   void simulate3D(std::vector<double> &c);
   inline double getBCFromFlux(boundary_condition bc, double nearest_value,
                               double neighbor_alpha);
+    inline void solveLES();
   void updateInternals();
 
   std::vector<boundary_condition> bc;