diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp
index 2323531..d755e15 100644
--- a/examples/reference-FTCS_2D_closed.cpp
+++ b/examples/reference-FTCS_2D_closed.cpp
@@ -5,8 +5,8 @@
 using namespace std;
 
 int main(int argc, char *argv[]) {
-    int row = 11;
-    int col = 11;
+    int row = 50;
+    int col = 50;
     int domain_row = 10;
     int domain_col = 10;
 
@@ -45,14 +45,11 @@ int main(int argc, char *argv[]) {
     // Simulation
     Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
     sim.setTimestep(0.001);
-    sim.setIterations(7000);
-    sim.setOutputCSV(CSV_OUTPUT_ON);
-    sim.setOutputConsole(CONSOLE_OUTPUT_ON);
+    sim.setIterations(100);
+    sim.setOutputCSV(CSV_OUTPUT_OFF);
+    sim.setOutputConsole(CONSOLE_OUTPUT_OFF);
 
 
     // RUN
     sim.run();
-
-    cout << grid.getConcentrations() << endl;
-
 }
\ No newline at end of file
diff --git a/src/FTCS.cpp b/src/FTCS.cpp
index 4fc1498..fa9f546 100644
--- a/src/FTCS.cpp
+++ b/src/FTCS.cpp
@@ -9,6 +9,7 @@
 #include <cstddef>
 #include <tug/Boundary.hpp>
 #include <iostream>
+#include <omp.h>
 
 using namespace std;
 
@@ -277,6 +278,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 
     // inner cells
     // these are independent of the boundary condition type
+    omp_set_num_threads(10);
+    #pragma omp parallel for
     for (int row = 1; row < rowMax-1; row++) {
         for (int col = 1; col < colMax-1; col++) {
             concentrations_t1(row, col) = grid.getConcentrations()(row, col) 
@@ -296,6 +299,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
     // left without corners / looping over rows
     // hold column constant at index 0
     int col = 0;
+    #pragma omp parallel for
     for (int row = 1; row < rowMax-1; row++) {
         concentrations_t1(row, col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
@@ -312,6 +316,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
     // right without corners / looping over rows
     // hold column constant at max index
     col = colMax-1;
+    #pragma omp parallel for
     for (int row = 1; row < rowMax-1; row++) {
         concentrations_t1(row,col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
@@ -329,6 +334,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
     // top without corners / looping over columns
     // hold row constant at index 0
     int row = 0;
+    #pragma omp parallel for
     for (int col=1; col<colMax-1;col++){
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
             + timestep / (deltaRow*deltaRow) 
@@ -345,6 +351,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
     // bottom without corners / looping over columns
     // hold row constant at max index
     row = rowMax-1;
+    #pragma omp parallel for
     for(int col=1; col<colMax-1;col++){
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
             + timestep / (deltaRow*deltaRow) 
diff --git a/src/Simulation.cpp b/src/Simulation.cpp
index d07065f..5fd552c 100644
--- a/src/Simulation.cpp
+++ b/src/Simulation.cpp
@@ -154,7 +154,7 @@ void Simulation::run() {
     }
 
     if (approach == FTCS_APPROACH) {
-        
+        auto begin = std::chrono::high_resolution_clock::now();
         for (int i = 0; i < iterations; i++) {
             if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
                 printConcentrationsConsole();
@@ -165,6 +165,9 @@ void Simulation::run() {
 
             FTCS(grid, bc, timestep);
         }
+        auto end = std::chrono::high_resolution_clock::now();
+        auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
+        std::cout << milliseconds.count() << endl;
 
     } else if (approach == BTCS_APPROACH) {