diff --git a/src/diffusion.cpp b/src/diffusion.cpp
index e668388..11410ff 100644
--- a/src/diffusion.cpp
+++ b/src/diffusion.cpp
@@ -2,6 +2,13 @@
 
 #include <Eigen/SparseLU>
 #include <Eigen/src/Core/Matrix.h>
+#include <Eigen/src/Core/util/Constants.h>
+#include <Eigen/src/OrderingMethods/Ordering.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 <iostream>
 
 void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
             double timestep) {
@@ -11,7 +18,8 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
 
   int size = x * y;
 
-  Eigen::VectorXd d = Eigen::VectorXd::Constant(size, 0);
+  Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
+  Eigen::VectorXd x_out(size);
   std::vector<T> tripletList;
   tripletList.reserve(((x - 1) * (y - 1) * 5));
 
@@ -25,6 +33,21 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
       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));
+
+      b[(i * x) + j] = -c[(i * x) + j];
     }
   }
+
+  Eigen::SparseMatrix<double> A(size, size);
+  A.setFromTriplets(tripletList.begin(), tripletList.end());
+
+  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
+      solver;
+  solver.analyzePattern(A);
+
+  solver.factorize(A);
+
+  std::cout << A << std::endl;
+
+  x_out = solver.solve(b);
 }
diff --git a/src/diffusion.hpp b/src/diffusion.hpp
index 5e04332..6645f8c 100644
--- a/src/diffusion.hpp
+++ b/src/diffusion.hpp
@@ -6,6 +6,6 @@
 
 typedef Eigen::Triplet<double> T;
 
-extern void BTCS2D(int x, int y, std::vector<double> &c, double alpha,
-                   double timestep);
+extern void BTCS2D(int x, int y, std::vector<double> &c,
+                   std::vector<double> &alpha, double timestep);
 #endif // DIFFUSION_H_