diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 9094a2c..1c38bf6 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -5,6 +5,8 @@
 
 #include <Eigen/src/Core/Map.h>
 #include <Eigen/src/Core/Matrix.h>
+#include <Eigen/src/SparseCore/SparseMatrix.h>
+#include <Eigen/src/SparseCore/SparseMatrixBase.h>
 #include <algorithm>
 #include <array>
 #include <cassert>
@@ -86,28 +88,32 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
                                              int size,
                                              const DVectorRowMajor &t0_c) {
 
-  reserveMemory(size, BTCS_MAX_DEP_PER_CELL);
+  Eigen::SparseMatrix<double> A_matrix;
+  Eigen::VectorXd b_vector;
+  Eigen::VectorXd x_vector;
 
-  fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step);
-  fillVectorFromRow(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, dx,
-                    time_step);
+  A_matrix.resize(size + 2, size + 2);
+  A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL));
 
-  solveLES();
+  b_vector.resize(size + 2);
+  x_vector.resize(size + 2);
+
+  fillMatrixFromRow(A_matrix, alpha.row(0), bc.row(0), size, dx, time_step);
+  fillVectorFromRow(b_vector, c, alpha, bc, Eigen::VectorXd::Constant(size, 0),
+                    size, dx, time_step);
+
+  // 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);
 
   c = x_vector.segment(1, size);
 }
 
-inline void Diffusion::BTCSDiffusion::reserveMemory(int size,
-                                                    int max_count_per_line) {
-  size += 2;
-
-  A_matrix.resize(size, size);
-  A_matrix.reserve(Eigen::VectorXi::Constant(size, max_count_per_line));
-
-  b_vector.resize(size);
-  x_vector.resize(size);
-}
-
 void Diffusion::BTCSDiffusion::simulate1D(
     Eigen::Map<DVectorRowMajor> &c, Eigen::Map<const DVectorRowMajor> &alpha,
     Eigen::Map<const BCVectorRowMajor> &bc) {
@@ -208,9 +214,9 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
   return t0_c;
 }
 
-inline void Diffusion::BTCSDiffusion::fillMatrixFromRow(
-    const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size,
-    double dx, double time_step) {
+void Diffusion::BTCSDiffusion::fillMatrixFromRow(
+    Eigen::SparseMatrix<double> &A_matrix, const DVectorRowMajor &alpha,
+    const BCVectorRowMajor &bc, int size, double dx, double time_step) {
 
   Diffusion::boundary_condition left = bc[0];
   Diffusion::boundary_condition right = bc[size - 1];
@@ -247,10 +253,10 @@ inline void Diffusion::BTCSDiffusion::fillMatrixFromRow(
   }
 }
 
-inline void Diffusion::BTCSDiffusion::fillVectorFromRow(
-    const DVectorRowMajor &c, const DVectorRowMajor &alpha,
-    const BCVectorRowMajor &bc, const DVectorRowMajor &t0_c, int size,
-    double dx, double time_step) {
+void Diffusion::BTCSDiffusion::fillVectorFromRow(
+    Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
+    const DVectorRowMajor &alpha, const BCVectorRowMajor &bc,
+    const DVectorRowMajor &t0_c, int size, double dx, double time_step) {
 
   Diffusion::boundary_condition left = bc[0];
   Diffusion::boundary_condition right = bc[size - 1];
@@ -341,14 +347,3 @@ inline auto Diffusion::BTCSDiffusion::getBCFromFlux(boundary_condition bc,
 
   return val;
 }
-
-inline void Diffusion::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 f479068..ef4c4d6 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -7,6 +7,7 @@
 #include <Eigen/src/Core/Map.h>
 #include <Eigen/src/Core/Matrix.h>
 #include <Eigen/src/Core/util/Constants.h>
+#include <Eigen/src/SparseCore/SparseMatrix.h>
 #include <cstddef>
 #include <tuple>
 #include <type_traits>
@@ -134,28 +135,24 @@ private:
                  const BCMatrixRowMajor &bc, double time_step, double dx)
       -> DMatrixRowMajor;
 
-  inline void fillMatrixFromRow(const DVectorRowMajor &alpha,
-                                const BCVectorRowMajor &bc, int size, double dx,
-                                double time_step);
-  inline void fillVectorFromRow(const DVectorRowMajor &c,
-                                const DVectorRowMajor &alpha,
-                                const BCVectorRowMajor &bc,
-                                const DVectorRowMajor &t0_c, int size,
-                                double dx, double time_step);
+  void fillMatrixFromRow(Eigen::SparseMatrix<double> &A_matrix,
+                         const DVectorRowMajor &alpha,
+                         const BCVectorRowMajor &bc, int size, double dx,
+                         double time_step);
+
+  void fillVectorFromRow(Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
+                         const DVectorRowMajor &alpha,
+                         const BCVectorRowMajor &bc,
+                         const DVectorRowMajor &t0_c, int size, double dx,
+                         double time_step);
   void simulate3D(std::vector<double> &c);
 
-  inline void reserveMemory(int size, int max_count_per_line);
   inline static auto getBCFromFlux(Diffusion::boundary_condition bc,
                                    double neighbor_c, double neighbor_alpha)
       -> double;
 
-  void solveLES();
   void updateInternals();
 
-  Eigen::SparseMatrix<double> A_matrix;
-  Eigen::VectorXd b_vector;
-  Eigen::VectorXd x_vector;
-
   double time_step;
   unsigned int grid_dim;