diff --git a/CMakeLists.txt b/CMakeLists.txt
index 068023f..0fe7fb6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,6 +6,9 @@ project(Diffusion CXX)
 set(CMAKE_CXX_STANDARD 14)
 
 find_package(Eigen3 REQUIRED NO_MODULE)
+find_package(OpenMP)
+
+option(USE_OPENMP "Compile with OpenMP support" ON)
 
 add_subdirectory(app)
 add_subdirectory(src)
diff --git a/app/Rcpp-BTCS-1d.cpp b/app/Rcpp-BTCS-1d.cpp
new file mode 100644
index 0000000..8949d85
--- /dev/null
+++ b/app/Rcpp-BTCS-1d.cpp
@@ -0,0 +1,59 @@
+#include "../include/diffusion/BTCSDiffusion.hpp"
+#include "../include/diffusion/BoundaryCondition.hpp"
+#include <algorithm> // for copy, max
+#include <cmath>
+#include <iomanip>
+#include <iostream> // for std
+#include <vector>   // for vector
+#include <Rcpp.h>
+
+using namespace std;
+using namespace Diffusion;
+using namespace Rcpp;
+//using namespace Eigen;
+
+// [[Rcpp::depends(RcppEigen)]]
+// [[Rcpp::plugins("cpp11")]]
+// [[Rcpp::export]]
+std::vector<double> & diff1D(int n,
+                             double length,
+                             std::vector<double> & field,
+                             std::vector<double> & alpha,
+                             double timestep,
+                             double bc_left,
+                             double bc_right,
+                             int iterations) {
+    // dimension of grid
+    int dim = 1;
+    
+    // create input + diffusion coefficients for each grid cell
+    // std::vector<double> alpha(n, 1 * pow(10, -1));
+    // std::vector<double> field(n, 1 * std::pow(10, -6));
+    std::vector<boundary_condition> bc(n, {0,0});
+    
+    // create instance of diffusion module
+    BTCSDiffusion diffu(dim);
+    
+    diffu.setXDimensions(length, n);
+    
+    // set the boundary condition for the left ghost cell to dirichlet
+    bc[0] = {Diffusion::BC_CONSTANT, bc_left};
+    bc[n-1] = {Diffusion::BC_CONSTANT, bc_right};
+    
+    // set timestep for simulation to 1 second
+    diffu.setTimestep(timestep);
+    
+    //cout << setprecision(12);
+    
+    // loop 100 times
+    // output is currently generated by the method itself
+    for (int i = 0; i < iterations; i++) {
+        diffu.simulate(field.data(), alpha.data(), bc.data());
+    }
+
+    // for (auto & cell : field) {
+    //     Rcout << cell << "\n";
+    // }
+    
+    return(field);
+}
diff --git a/app/Rcpp-BTCS-2d.cpp b/app/Rcpp-BTCS-2d.cpp
new file mode 100644
index 0000000..3a5e404
--- /dev/null
+++ b/app/Rcpp-BTCS-2d.cpp
@@ -0,0 +1,60 @@
+#include "../include/diffusion/BTCSDiffusion.hpp"
+#include "../include/diffusion/BoundaryCondition.hpp"
+#include <algorithm> // for copy, max
+#include <cmath>
+#include <iomanip>
+#include <iostream> // for std
+#include <vector>   // for vector
+#include <Rcpp.h>
+
+using namespace std;
+using namespace Diffusion;
+using namespace Rcpp;
+//using namespace Eigen;
+
+// [[Rcpp::depends(RcppEigen)]]
+// [[Rcpp::plugins("cpp11")]]
+// [[Rcpp::export]]
+std::vector<double> & diff2D(int nx,
+                             int ny,
+                             double lenx,
+                             double leny,
+                             std::vector<double> & field,
+                             std::vector<double> & alpha,
+                             double timestep,
+                             int iterations)
+{
+    // problem dimensionality
+    int dim = 2;
+    // total number of grid cells
+    int n = nx*ny;
+    
+    std::vector<boundary_condition> bc(n, {0,0});
+    
+    // create instance of diffusion module
+    BTCSDiffusion diffu(dim);
+    
+    diffu.setXDimensions(lenx, nx);
+    diffu.setXDimensions(leny, ny);
+    
+    // set the boundary condition for the left ghost cell to dirichlet
+    // bc[0] = {Diffusion::BC_CONSTANT, bc_left};
+    // bc[n-1] = {Diffusion::BC_CONSTANT, bc_right};
+    
+    // set timestep for simulation to 1 second
+    diffu.setTimestep(timestep);
+    
+    //cout << setprecision(12);
+    
+    // loop 100 times
+    // output is currently generated by the method itself
+    for (int i = 0; i < iterations; i++) {
+        diffu.simulate(field.data(), alpha.data(), bc.data());
+    }
+
+    // for (auto & cell : field) {
+    //     Rcout << cell << "\n";
+    // }
+    
+    return(field);
+}
diff --git a/app/Rcpp-interface.R b/app/Rcpp-interface.R
new file mode 100644
index 0000000..fc1280b
--- /dev/null
+++ b/app/Rcpp-interface.R
@@ -0,0 +1,159 @@
+## Time-stamp: "Last modified 2022-03-16 14:01:11 delucia"
+library(Rcpp)
+library(RcppEigen)
+library(ReacTran)
+library(deSolve)
+
+options(width=110)
+
+setwd("app")
+
+## This creates the "diff1D" function with our BTCSdiffusion code
+sourceCpp("Rcpp-BTCS-1d.cpp")
+
+### FTCS explicit (same name)
+sourceCpp("RcppFTCS.cpp")
+
+
+
+## Grid 101
+## Set initial conditions
+N <- 1001
+D.coeff <- 1E-3
+C0 <- 1             ## Initial concentration (mg/L)
+X0 <- 0             ## Location of initial concentration (m)
+## Yini <- c(C0, rep(0,N-1))
+
+## Ode1d solution
+xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N)
+x <- xgrid$x.mid
+Diffusion <- function (t, Y, parms){
+  tran <- tran.1D(C = Y, C.up = 0, C.down = 0, D = parms$D, dx = xgrid)
+  return(list(tran$dC))
+}
+
+
+## gaussian pulse as initial condition
+sigma <- 0.02
+Yini <- 0.5*exp(-0.5*((x-1/2.0)**2)/sigma**2)
+
+## plot(x, Yini, type="l")
+
+parms1 <- list(D=D.coeff)
+# 1 timestep, 10 s
+times <- seq(from = 0, to = 1, by = 0.1)
+
+system.time({
+    out1 <- ode.1D(y = Yini, times = times, func = Diffusion,
+                  parms = parms1, dimens = N)[11,-1]
+})
+
+## Now with BTCS 
+alpha <- rep(D.coeff, N)
+
+system.time({
+    out2 <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = 0.1, 0, 0, iterations = 10)
+})
+
+## plot(out1, out2)
+## abline(0,1)
+
+## matplot(cbind(out1,out2),type="l", col=c("black","red"),lty="solid", lwd=2,
+##         xlab="grid element", ylab="Concentration", las=1)
+## legend("topright", c("ReacTran ode1D", "BTCS 1d"), text.col=c("black","red"), bty = "n")
+
+
+
+
+
+system.time({
+    out3 <- RcppFTCS(n=N, length=1, field=Yini, alpha=1E-3, bc_left = 0, bc_right = 0, timestep = 1)
+})
+
+## Poor man's 
+mm <- colMeans(rbind(out2,out3))
+
+
+
+matplot(cbind(Yini,out1, out2, out3, mm),type="l", col=c("grey","black","red","blue","green4"), lty="solid", lwd=2,
+        xlab="grid element", ylab="Concentration", las=1)
+legend("topright", c("init","ReacTran ode1D", "BTCS 1d", "FTCS", "poor man's CN"), text.col=c("grey","black","red","blue","green4"), bty = "n")
+
+
+sum(Yini)
+sum(out1)
+sum(out2)
+sum(out3)
+sum(mm)
+
+## Yini <- 0.2*sin(pi/0.1*x)+0.2
+## plot(Yini)
+
+## plot(out3)
+
+Fun <- function(dx) {
+    tmp <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = floor(1/dx))
+    sqrt(sum({out1-tmp}^2))
+}
+
+reso <- optimise(f=Fun, interval=c(1E-5, 1E-1), maximum = FALSE)
+
+
+dx <- 0.0006038284
+floor(1/dx)
+
+1/dx
+
+system.time({
+    out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1656)
+})
+
+matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2,
+        xlab="grid element", ylab="Concentration", las=1)
+legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n")
+
+
+dx <- 0.05
+
+system.time({
+    out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1/dx)
+})
+
+matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2,
+        xlab="grid element", ylab="Concentration", las=1)
+legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n")
+
+Matplot
+
+
+
+
+## This creates the "diff1D" function with our BTCSdiffusion code
+sourceCpp("Rcpp-BTCS-2d.cpp")
+
+n <- 256
+a2d <- rep(1E-3, n^2)
+
+init2d <- readRDS("gs1.rds")
+
+ll <- {init2d - min(init2d)}/diff(range(init2d))
+
+system.time({
+    res1 <- diff2D(nx=N, ny=N, lenx=1, leny=1, field=ll, alpha=a2d, timestep = 0.1, iterations = 10)
+})
+
+hist(ll,32)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/app/Rcpp-interface.cpp b/app/Rcpp-interface.cpp
new file mode 100644
index 0000000..8949d85
--- /dev/null
+++ b/app/Rcpp-interface.cpp
@@ -0,0 +1,59 @@
+#include "../include/diffusion/BTCSDiffusion.hpp"
+#include "../include/diffusion/BoundaryCondition.hpp"
+#include <algorithm> // for copy, max
+#include <cmath>
+#include <iomanip>
+#include <iostream> // for std
+#include <vector>   // for vector
+#include <Rcpp.h>
+
+using namespace std;
+using namespace Diffusion;
+using namespace Rcpp;
+//using namespace Eigen;
+
+// [[Rcpp::depends(RcppEigen)]]
+// [[Rcpp::plugins("cpp11")]]
+// [[Rcpp::export]]
+std::vector<double> & diff1D(int n,
+                             double length,
+                             std::vector<double> & field,
+                             std::vector<double> & alpha,
+                             double timestep,
+                             double bc_left,
+                             double bc_right,
+                             int iterations) {
+    // dimension of grid
+    int dim = 1;
+    
+    // create input + diffusion coefficients for each grid cell
+    // std::vector<double> alpha(n, 1 * pow(10, -1));
+    // std::vector<double> field(n, 1 * std::pow(10, -6));
+    std::vector<boundary_condition> bc(n, {0,0});
+    
+    // create instance of diffusion module
+    BTCSDiffusion diffu(dim);
+    
+    diffu.setXDimensions(length, n);
+    
+    // set the boundary condition for the left ghost cell to dirichlet
+    bc[0] = {Diffusion::BC_CONSTANT, bc_left};
+    bc[n-1] = {Diffusion::BC_CONSTANT, bc_right};
+    
+    // set timestep for simulation to 1 second
+    diffu.setTimestep(timestep);
+    
+    //cout << setprecision(12);
+    
+    // loop 100 times
+    // output is currently generated by the method itself
+    for (int i = 0; i < iterations; i++) {
+        diffu.simulate(field.data(), alpha.data(), bc.data());
+    }
+
+    // for (auto & cell : field) {
+    //     Rcout << cell << "\n";
+    // }
+    
+    return(field);
+}
diff --git a/app/RcppFTCS.cpp b/app/RcppFTCS.cpp
new file mode 100644
index 0000000..dbe473b
--- /dev/null
+++ b/app/RcppFTCS.cpp
@@ -0,0 +1,40 @@
+//  Time-stamp: "Last modified 2022-03-15 18:15:39 delucia"
+#include <Rcpp.h>
+#include <iostream> // for std
+#include <vector>   // for vector
+
+using namespace std;
+using namespace Rcpp;
+
+// [[Rcpp::plugins("cpp11")]]
+// [[Rcpp::export]]
+NumericVector RcppFTCS(int n,
+                       double length,    
+                       NumericVector & field,
+                       double alpha,
+                       double bc_left,
+                       double bc_right,
+                       double timestep)
+{
+    // dimension of grid
+    NumericVector ext (clone(field));
+    double dx = length / ((double) n - 1.);
+    double dt = 0.25*dx*dx/alpha;
+
+    
+    double afac = alpha*dt/dx/dx;
+    int iter = (int) (timestep/dt);
+
+    Rcout << "dt: " << dt << "; inner iterations: " << iter << endl;
+
+
+    for (int it = 0; it < iter; it++){
+        for (int i = 1; i < ext.size()-1; i++) {
+            ext[i] = (1. - 2*afac)*ext[i] + afac*(ext[i+1]+ext[i-1]);
+        }
+        ext[0] = bc_left;
+        ext[n-1] = bc_right;
+    }
+    return(ext);
+}
+
diff --git a/app/main_1D.cpp b/app/main_1D.cpp
index f7ec59b..b2423f2 100644
--- a/app/main_1D.cpp
+++ b/app/main_1D.cpp
@@ -1,5 +1,6 @@
-#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
-#include "BoundaryCondition.hpp"
+#include <diffusion/BTCSDiffusion.hpp>
+#include <diffusion/BoundaryCondition.hpp>
+
 #include <algorithm> // for copy, max
 #include <cmath>
 #include <iomanip>
diff --git a/app/main_2D.cpp b/app/main_2D.cpp
index e67004c..6a35bba 100644
--- a/app/main_2D.cpp
+++ b/app/main_2D.cpp
@@ -1,5 +1,5 @@
-#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
-#include "BoundaryCondition.hpp"
+#include <diffusion/BTCSDiffusion.hpp>
+#include <diffusion/BoundaryCondition.hpp>
 #include <algorithm>         // for copy, max
 #include <cmath>
 #include <iomanip>
diff --git a/app/main_2D_mdl.cpp b/app/main_2D_mdl.cpp
index 900741b..c040cfe 100644
--- a/app/main_2D_mdl.cpp
+++ b/app/main_2D_mdl.cpp
@@ -1,5 +1,5 @@
-#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
-#include "BoundaryCondition.hpp"
+#include <diffusion/BTCSDiffusion.hpp>
+#include <diffusion/BoundaryCondition.hpp>
 #include <algorithm>         // for copy, max
 #include <cmath>
 #include <iomanip>
diff --git a/src/BTCSDiffusion.hpp b/include/diffusion/BTCSDiffusion.hpp
similarity index 86%
rename from src/BTCSDiffusion.hpp
rename to include/diffusion/BTCSDiffusion.hpp
index f479068..ef4c4d6 100644
--- a/src/BTCSDiffusion.hpp
+++ b/include/diffusion/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;
 
diff --git a/src/BoundaryCondition.hpp b/include/diffusion/BoundaryCondition.hpp
similarity index 100%
rename from src/BoundaryCondition.hpp
rename to include/diffusion/BoundaryCondition.hpp
diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 9094a2c..8fdd020 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -1,10 +1,12 @@
-#include "BTCSDiffusion.hpp"
-#include "BoundaryCondition.hpp"
+#include "diffusion/BTCSDiffusion.hpp"
+#include "diffusion/BoundaryCondition.hpp"
 
 #include <Eigen/SparseLU>
 
 #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>
@@ -16,6 +18,12 @@
 #include <tuple>
 #include <vector>
 
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_thread_num() 0
+#endif
+
 #include <iostream>
 
 constexpr int BTCS_MAX_DEP_PER_CELL = 3;
@@ -86,28 +94,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) {
@@ -137,6 +149,7 @@ void Diffusion::BTCSDiffusion::simulate2D(
 
   t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
 
+#pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_rows; i++) {
     DVectorRowMajor input_field = c.row(i);
     simulate_base(input_field, bc.row(i), alpha.row(i), dx, local_dt, n_cols,
@@ -149,6 +162,7 @@ void Diffusion::BTCSDiffusion::simulate2D(
   t0_c =
       calc_t0_c(c.transpose(), alpha.transpose(), bc.transpose(), local_dt, dx);
 
+#pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_cols; i++) {
     DVectorRowMajor input_field = c.col(i);
     simulate_base(input_field, bc.col(i), alpha.col(i), dx, local_dt, n_rows,
@@ -180,7 +194,8 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
                  (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
   }
 
-  // then iterate over inlet
+// then iterate over inlet
+#pragma omp parallel for private(y_values) schedule(dynamic)
   for (int i = 1; i < n_rows - 1; i++) {
     for (int j = 0; j < n_cols; j++) {
 
@@ -208,9 +223,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 +262,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 +356,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/CMakeLists.txt b/src/CMakeLists.txt
index ffc3958..30a0f62 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,12 @@
-add_library(diffusion OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp)
+set(HEADER_LIST "${Diffusion_SOURCE_DIR}/include/diffusion/BTCSDiffusion.hpp"
+  "${Diffusion_SOURCE_DIR}/include/diffusion/BoundaryCondition.hpp")
+
+add_library(diffusion STATIC BTCSDiffusion.cpp ${HEADER_LIST})
+
 target_link_libraries(diffusion Eigen3::Eigen)
-target_include_directories(diffusion PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+
+if(USE_OPENMP AND OpenMP_CXX_FOUND)
+  target_link_libraries(diffusion OpenMP::OpenMP_CXX)
+endif()
+
+target_include_directories(diffusion PUBLIC ../include)