diff --git a/src/diffusion.cpp b/src/diffusion.cpp
index 8d9ca5f..7ecad7f 100644
--- a/src/diffusion.cpp
+++ b/src/diffusion.cpp
@@ -4,10 +4,14 @@
 #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/SparseCholesky>
+#include<Eigen/SparseQR>
+#include <Eigen/src/SparseQR/SparseQR.h>
 #include <iostream>
 
 void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
@@ -16,7 +20,7 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
   double dx = 1. / x;
   double dy = 1. / y;
 
-  int size = (x - 2) * (y - 2);
+  int size = (x * y) - 4;
 
   int local_x = x - 2;
 
@@ -25,22 +29,62 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
   std::vector<T> tripletList;
   tripletList.reserve(size * 5);
 
-  for (int i = 0; i < (size-local_x); i++) {
-    int current = local_x + i;
-    double sx = (alpha[current] * timestep) / (dx * dx);
-    double sy = (alpha[current] * timestep) / (dy * dy);
+  for (int i = x - 1 ; i < 2*x - 3 ; i++) {
+    double sx = (alpha[i + 2] * timestep) / (dx * dx);
+    double sy = (alpha[i + 2] * timestep) / (dy * dy);
 
-    tripletList.push_back(T(i, current, (1 + 2 * sx + 2 * sy)));
-    tripletList.push_back(T(i, current + local_x, sy));
-    tripletList.push_back(T(i, current - local_x, sy));
-    tripletList.push_back(T(i, current + 1, sx));
-    tripletList.push_back(T(i, current - 1, sx));
 
-    std::cout << current << std::endl;
+    tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
+    tripletList.push_back(T(i, i - (x - 1), sy));
+    tripletList.push_back(T(i, i + x, sy));
+    tripletList.push_back(T(i, i + 1, sx));
+    tripletList.push_back(T(i, i - 1, sx));
 
-    b[i] = -c[x+i];
+    b[i] = -c[i+2];
   }
 
+  for (int i = 2*x - 1; i < (y-2)*x-3; i++) {
+    double sx = (alpha[i + 2] * timestep) / (dx * dx);
+    double sy = (alpha[i + 2] * timestep) / (dy * dy);
+
+    tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
+    tripletList.push_back(T(i, i - x, sy));
+    tripletList.push_back(T(i, i + x, sy));
+    tripletList.push_back(T(i, i + 1, sx));
+    tripletList.push_back(T(i, i - 1, sx));
+
+    b[i] = -c[i+2];
+  }
+
+  for (int i = (y-2)*x-1; i < (y-1)*x-3; i++) {
+    double sx = (alpha[i + 2] * timestep) / (dx * dx);
+    double sy = (alpha[i + 2] * timestep) / (dy * dy);
+
+    tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
+    tripletList.push_back(T(i, i - x, sy));
+    tripletList.push_back(T(i, i + (x-1), sy));
+    tripletList.push_back(T(i, i + 1, sx));
+    tripletList.push_back(T(i, i - 1, sx));
+
+    b[i] = -c[i+2];
+  }
+
+  // for (int i = 0; i < (size-local_x); i++) {
+  //   int current = local_x + i;
+  //   double sx = (alpha[current] * timestep) / (dx * dx);
+  //   double sy = (alpha[current] * timestep) / (dy * dy);
+
+  //   tripletList.push_back(T(i, current, (1 + 2 * sx + 2 * sy)));
+  //   tripletList.push_back(T(i, current + local_x, sy));
+  //   tripletList.push_back(T(i, current - local_x, sy));
+  //   tripletList.push_back(T(i, current + 1, sx));
+  //   tripletList.push_back(T(i, current - 1, sx));
+
+  //   std::cout << current << std::endl;
+
+  //   b[i] = -c[x+i];
+  // }
+
   // for (int i = 0; i < y; i++) {
   //   for (int j = 0; j < x; j++) {
   //     double sx = (alpha[i * y + j] * timestep) / (dx * dx);
@@ -56,16 +100,23 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
   //   }
   // }
 
-  Eigen::SparseMatrix<double> A(size, (x*y)-4);
+  std::cout << b << std::endl;
+
+  Eigen::SparseMatrix<double> A(size, (x * y) - 4);
   A.setFromTriplets(tripletList.begin(), tripletList.end());
 
-  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
-      solver;
+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 << A << std::endl;
+  std::cout << solver.lastErrorMessage() << std::endl;
 
   x_out = solver.solve(b);
+
+  std::cout << x_out << std::endl;
 }
diff --git a/src/main.cpp b/src/main.cpp
index 6c18b6f..8269d9b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -7,12 +7,14 @@ using namespace std;
 
 int main(int argc, char *argv[]) {
 
-  int x = 10;
-  int y = 10;
+  int x = 5;
+  int y = 5;
 
-  std::vector<double> alpha(x * y, 1.5 * pow(10, -9));
+  std::vector<double> alpha(x * y, -1.5 * pow(10, -2));
   std::vector<double> input(x * y, 0);
-  input[x + 1] = 5 * std::pow(10, 6);
+  input[x + 1] = 5.55 * std::pow(10, -6);
+  input[x + 2] = 5.5556554 * std::pow(10, -6);
+  input[x + 3] = 5.234564213 * std::pow(10, -6);
 
   BTCS2D(x, y, input, alpha, 10.);