diff --git a/src/diffusion.cpp b/src/diffusion.cpp
new file mode 100644
index 0000000..e668388
--- /dev/null
+++ b/src/diffusion.cpp
@@ -0,0 +1,30 @@
+#include "diffusion.hpp"
+
+#include <Eigen/SparseLU>
+#include <Eigen/src/Core/Matrix.h>
+
+void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
+            double timestep) {
+
+  double dx = 1. / x;
+  double dy = 1. / y;
+
+  int size = x * y;
+
+  Eigen::VectorXd d = Eigen::VectorXd::Constant(size, 0);
+  std::vector<T> tripletList;
+  tripletList.reserve(((x - 1) * (y - 1) * 5));
+
+  for (int i = 1; i < y - 1; i++) {
+    for (int j = 1; j < x - 1; j++) {
+      double sx = (alpha[i * y + j] * timestep) / (dx * dx);
+      double sy = (alpha[i * y + j] * timestep) / (dy * dy);
+
+      tripletList.push_back(T((i * x) + j, (i * x) + j, (1 + 2 * sx + 2 * sy)));
+      tripletList.push_back(T((i * x) + j, ((i * x) + j) + x, sy));
+      tripletList.push_back(T((i * x) + j, ((i * x) + j) - x, sy));
+      tripletList.push_back(T((i * x) + j, ((i * x) + j) + 1, sx));
+      tripletList.push_back(T((i * x) + j, ((i * x) + j) - 1, sx));
+    }
+  }
+}
diff --git a/src/diffusion.hpp b/src/diffusion.hpp
new file mode 100644
index 0000000..5e04332
--- /dev/null
+++ b/src/diffusion.hpp
@@ -0,0 +1,11 @@
+#ifndef DIFFUSION_H_
+#define DIFFUSION_H_
+
+#include <Eigen/SparseCore>
+#include <vector>
+
+typedef Eigen::Triplet<double> T;
+
+extern void BTCS2D(int x, int y, std::vector<double> &c, double alpha,
+                   double timestep);
+#endif // DIFFUSION_H_