diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
new file mode 100644
index 0000000..919a1e5
--- /dev/null
+++ b/src/BTCSDiffusion.cpp
@@ -0,0 +1,96 @@
+#include "BTCSDiffusion.hpp"
+
+#include <Eigen/SparseCholesky>
+#include <Eigen/SparseLU>
+#include <Eigen/SparseQR>
+#include <Eigen/src/Core/Matrix.h>
+#include <Eigen/src/Core/util/Constants.h>
+#include <Eigen/src/OrderingMethods/Ordering.h>
+#include <Eigen/src/SparseCholesky/SimplicialCholesky.h>
+#include <Eigen/src/SparseCore/SparseMap.h>
+#include <Eigen/src/SparseCore/SparseMatrix.h>
+#include <Eigen/src/SparseCore/SparseMatrixBase.h>
+#include <Eigen/src/SparseLU/SparseLU.h>
+#include <Eigen/src/SparseQR/SparseQR.h>
+
+#include <iomanip>
+#include <iostream>
+
+const BCSide BTCSDiffusion::LEFT = 0;
+const BCSide BTCSDiffusion::RIGHT = 1;
+
+BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { this->grid_dim = 1; }
+BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) {
+
+  this->grid_dim = 2;
+}
+BTCSDiffusion::BTCSDiffusion(int x, int y, int z)
+    : dim_x(x), dim_y(y), dim_z(z) {
+
+  this->grid_dim = 3;
+}
+
+void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
+                             double timestep) {
+  double dx = 1. / this->dim_x;
+
+  int size = this->dim_x + 2;
+
+  Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
+  Eigen::VectorXd x_out(size);
+  std::vector<T> tripletList;
+  tripletList.reserve(c.size() * 3 + bc.size());
+
+  int A_line = 0;
+
+  for (int i = 1; i < this->dim_x + 1; i++) {
+    double sx = (alpha[i - 1] * timestep) / (dx * dx);
+
+    tripletList.push_back(T(A_line, i, (-1. - 2. * sx)));
+
+    tripletList.push_back(T(A_line, i - 1, sx));
+    tripletList.push_back(T(A_line, i + 1, sx));
+
+    b[A_line] = -c[i - 1];
+    A_line++;
+  }
+
+  tripletList.push_back(T(A_line, 0, 1));
+  if (bc[0] == -1)
+    b[A_line] = c[0];
+  else
+    b[A_line] = this->bc[0];
+
+  A_line++;
+  tripletList.push_back(T(A_line, size - 1, 1));
+  // b[A_line] = bc[1];
+  if (bc[0] == -1)
+    b[A_line] = c[c.size() - 1];
+  else
+    b[A_line] = this->bc[1];
+
+  // std::cout << b << std::endl;
+
+  Eigen::SparseMatrix<double> A(size, size);
+  A.setFromTriplets(tripletList.begin(), tripletList.end());
+
+  // std::cout << A << std::endl;
+  Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
+      solver;
+
+  // Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
+  //     solver;
+  solver.analyzePattern(A);
+
+  solver.factorize(A);
+
+  std::cout << solver.lastErrorMessage() << std::endl;
+
+  x_out = solver.solve(b);
+
+  std::cout << std::setprecision(10) << x_out << std::endl << std::endl;
+
+  for (int i = 0; i < c.size(); i++) {
+    c[i] = x_out[i + 1];
+  }
+}
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
new file mode 100644
index 0000000..cee17c9
--- /dev/null
+++ b/src/BTCSDiffusion.hpp
@@ -0,0 +1,32 @@
+#ifndef BTCSDIFFUSION_H_
+#define BTCSDIFFUSION_H_
+
+#include <vector>
+#include <Eigen/SparseCore>
+
+typedef int BCSide;
+typedef Eigen::Triplet<double> T;
+
+class BTCSDiffusion {
+
+public:
+  static const BCSide LEFT;
+  static const BCSide RIGHT;
+
+  BTCSDiffusion(int x);
+  BTCSDiffusion(int x, int y);
+  BTCSDiffusion(int x, int y, int z);
+
+  void setBoundaryCondition(std::vector<double> input, BCSide side);
+  void simulate(std::vector<double> &c, std::vector<double> &alpha,
+                double timestep);
+
+private:
+  std::vector<double> bc;
+  int grid_dim;
+  int dim_x;
+  int dim_y;
+  int dim_z;
+};
+
+#endif // BTCSDIFFUSION_H_
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 33e2d28..97ce3f1 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,5 +1,8 @@
 add_library(diffusion OBJECT diffusion.cpp diffusion.hpp)
 target_link_libraries(diffusion Eigen3::Eigen)
 
+add_library(diffusion_class OBJECT BTCSDiffusion.cpp)
+target_link_libraries(diffusion_class Eigen3::Eigen)
+
 add_executable(test main.cpp)
 target_link_libraries(test PUBLIC diffusion)