diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 1c38bf6..291c9d5 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -14,6 +14,7 @@
 #include <cstddef>
 #include <iomanip>
 #include <iterator>
+#include <omp.h>
 #include <ostream>
 #include <tuple>
 #include <vector>
@@ -143,6 +144,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,
@@ -155,6 +157,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,
@@ -186,7 +189,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++) {
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ffc3958..b12fd67 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,5 @@
+find_package(OpenMP)
+
 add_library(diffusion OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp)
-target_link_libraries(diffusion Eigen3::Eigen)
+target_link_libraries(diffusion Eigen3::Eigen OpenMP::OpenMP_CXX)
 target_include_directories(diffusion PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})