diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index d9f4baf..e6af65d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -15,6 +15,7 @@ build_release:
       - mkdir build && cd build
       - cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON ..
       - make -j$(nproc)
+      - cp ../test/FTCS_11_11_7000.csv test/
 
 test:
     stage: test
diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp
index 89c3b29..1deb307 100644
--- a/examples/reference-FTCS_2D_closed.cpp
+++ b/examples/reference-FTCS_2D_closed.cpp
@@ -1,4 +1,5 @@
 #include <tug/Simulation.hpp>
+#include "Eigen/Core"
 
 int main(int argc, char *argv[]) {
     int row = 11;
@@ -40,11 +41,13 @@ int main(int argc, char *argv[]) {
 
     // Simulation
     Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
-    //sim.setTimestep(0.001);
-    sim.setIterations(2);
-    sim.setOutputCSV(CSV_OUTPUT_VERBOSE);
+    sim.setTimestep(0.001);
+    sim.setIterations(7000);
+    sim.setOutputCSV(CSV_OUTPUT_ON);
+    sim.setOutputConsole(CONSOLE_OUTPUT_ON);
 
 
     // RUN
     sim.run();
+
 }
\ No newline at end of file
diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp
index d9658b0..91ef322 100644
--- a/include/tug/Simulation.hpp
+++ b/include/tug/Simulation.hpp
@@ -103,11 +103,11 @@ class Simulation {
         void printConcentrationsCSV(string filename);
 
         /**
-         * @brief Start the simulation with all of the previously set parameters. 
+         * @brief 
          * 
-         * @return auto 
+         * @return Grid 
          */
-        void run();
+        Grid run();
 
     private:
 
diff --git a/src/Simulation.cpp b/src/Simulation.cpp
index d74bc29..a975db3 100644
--- a/src/Simulation.cpp
+++ b/src/Simulation.cpp
@@ -25,12 +25,13 @@ Simulation::Simulation(Grid grid, Boundary bc, APPROACH approach) : grid(grid),
     double maxAlphaY = grid.getAlphaY().maxCoeff();
     double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
 
-    //double maxStableTimestep = minDelta / (2*maxAlpha); // Formula from Marco --> seems to be unstable
+    double maxStableTimestepMdl = minDelta / (2*maxAlpha); // Formula from Marco --> seems to be unstable
     double maxStableTimestep = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia
 
-    cout << maxStableTimestep << endl;
+    // cout << "Max stable time step MDL: " << maxStableTimestepMdl << endl;
+    // cout << "Max stable time step: " << maxStableTimestep << endl;
 
-    this->timestep = maxStableTimestep;
+    this->timestep = maxStableTimestep; 
 
     
     this->iterations = 1000;
@@ -145,7 +146,7 @@ void Simulation::printConcentrationsCSV(string filename) {
     file.close();
 }
 
-void Simulation::run() {
+Grid Simulation::run() {
     string filename;
     if (this->console_output > CONSOLE_OUTPUT_OFF) {
         printConcentrationsConsole();
@@ -189,4 +190,6 @@ void Simulation::run() {
     if (this->csv_output > CSV_OUTPUT_OFF) {
         printConcentrationsCSV(filename);
     }
+
+    return grid;
 }
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index ed46443..9431fc2 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -11,9 +11,16 @@ if(NOT DOCTEST_LIB)
   FetchContent_MakeAvailable(DocTest)
 endif()
 
-add_executable(testTug setup.cpp testBoundaryCondition.cpp testDiffusion.cpp)
+add_executable(testTug setup.cpp testBoundaryCondition.cpp testDiffusion.cpp testSimulation.cpp)
 target_link_libraries(testTug doctest tug)
 
+# get relative path of the CSV file
+get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH CACHE)
+# set relative path in header file
+configure_file(testSimulation.hpp.in testSimulation.hpp)
+# include test directory with generated header file from above
+target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}")
+
 add_custom_target(
   check
   COMMAND $<TARGET_FILE:testTug>
diff --git a/test/FTCS_11_11_7000.csv b/test/FTCS_11_11_7000.csv
new file mode 100644
index 0000000..7afe9bb
--- /dev/null
+++ b/test/FTCS_11_11_7000.csv
@@ -0,0 +1,13 @@
+0 0 0 0 0 0 0 0 0 0 0
+1.88664e-08 3.39962e-08 7.57021e-08 1.76412e-07 4.15752e-07 9.00973e-07 3.65403e-09 9.6579e-12 3.59442e-13 3.42591e-14 3.27595e-15
+1.19102e-06 1.95195e-06 3.92165e-06 8.30575e-06 1.78976e-05 3.60742e-05 2.02843e-07 2.24659e-09 2.35085e-10 2.64988e-11 2.90933e-12
+5.85009e-05 8.57948e-05 0.000151499 0.000284105 0.00054607 0.00100251 1.18494e-05 8.26706e-07 1.26394e-07 1.70309e-08 2.20525e-09
+0.00202345 0.00258511 0.00381783 0.00599829 0.00972689 0.0154873 0.0011152 0.000247309 5.10506e-05 8.64727e-06 1.37747e-06
+0.0205848 0.0217651 0.0238282 0.0262762 0.0285812 0.0303808 0.0374255 0.0204234 0.00674813 0.00160264 0.000338852
+0.0199112 0.0210265 0.0229587 0.0252019 0.0272007 0.0285225 0.0310896 0.0156681 0.00495399 0.00114176 0.000235573
+0.0184589 0.0194561 0.0211596 0.0230704 0.0246215 0.0253266 0.0228374 0.010038 0.00292986 0.000638799 0.000125942
+0.0166888 0.0175517 0.0190015 0.0205611 0.0216793 0.0218694 0.0160278 0.00593307 0.00155164 0.000312583 5.76964e-05
+0.0151262 0.0158758 0.0171155 0.0183949 0.019193 0.0190523 0.0115557 0.00356455 0.000817461 0.000148892 2.51893e-05
+0.0142177 0.0149034 0.0160255 0.0171522 0.0177843 0.0174891 0.0093846 0.0025221 0.000515469 8.51039e-05 1.31328e-05
+
+
diff --git a/test/TestUtils.cpp b/test/TestUtils.cpp
new file mode 100644
index 0000000..616b402
--- /dev/null
+++ b/test/TestUtils.cpp
@@ -0,0 +1,47 @@
+#include <ios>
+#include <iostream>
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <fstream>
+#include <sstream>
+#include <stdexcept>
+
+using namespace std;
+using namespace Eigen;
+
+MatrixXd CSV2Eigen(string file2Convert) {
+
+    vector<double> matrixEntries;
+
+    ifstream matrixDataFile(file2Convert);
+    if (matrixDataFile.fail()) {
+        throw invalid_argument("File probably non-existent!");
+    }
+
+    string matrixRowString;
+    string matrixEntry;
+    int matrixRowNumber = 0;
+
+    while(getline(matrixDataFile, matrixRowString)){
+        stringstream matrixRowStringStream(matrixRowString);
+        while(getline(matrixRowStringStream, matrixEntry, ' ')){
+            matrixEntries.push_back(stod(matrixEntry));
+        }
+        if (matrixRowString.length() > 1) {
+            matrixRowNumber++;
+        }
+    }
+
+    return Map<Matrix<double, Dynamic, Dynamic, RowMajor>>(matrixEntries.data(), matrixRowNumber, matrixEntries.size() / matrixRowNumber);
+}
+
+bool checkSimilarity(MatrixXd a, MatrixXd b, double precision=1e-5) {
+    return a.isApprox(b, precision);
+}
+
+bool checkSimilarityV2(MatrixXd a, MatrixXd b, double maxDiff) {
+
+    MatrixXd diff = a - b;
+    double maxCoeff = diff.maxCoeff();
+    return maxCoeff < maxDiff;
+}
\ No newline at end of file
diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp
new file mode 100644
index 0000000..8336cb7
--- /dev/null
+++ b/test/testSimulation.cpp
@@ -0,0 +1,67 @@
+#include <stdio.h>
+#include <doctest/doctest.h>
+#include <tug/Simulation.hpp>
+#include "TestUtils.cpp"
+#include <string>
+
+// include the configured header file
+#include <testSimulation.hpp>
+
+static Grid setupSimulation() {
+    int row = 11;
+    int col = 11;
+    int domain_row = 10;
+    int domain_col = 10;
+
+
+    // Grid
+    Grid grid = Grid(row, col);
+    grid.setDomain(domain_row, domain_col);
+
+    MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
+    concentrations(5,5) = 1;
+    grid.setConcentrations(concentrations);
+
+    MatrixXd alpha = MatrixXd::Constant(row, col, 1);
+    for (int i = 0; i < 5; i++) {
+        for (int j = 0; j < 6; j++) {
+            alpha(i, j) = 0.01;
+        }
+    }
+    for (int i = 0; i < 5; i++) {
+        for (int j = 6; j < 11; j++) {
+            alpha(i, j) = 0.001;
+        }
+    }
+    for (int i = 5; i < 11; i++) {
+        for (int j = 6; j < 11; j++) {
+            alpha(i, j) = 0.1;
+        }
+    }
+    grid.setAlpha(alpha, alpha);
+
+
+    // Boundary
+    Boundary bc = Boundary(grid);
+
+
+    // 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);
+
+
+    // RUN
+    return sim.run();
+}
+
+TEST_CASE("equality to reference matrix") {
+    // set string from the header file
+    string test_path = testSimulationCSVDir;
+    MatrixXd reference = CSV2Eigen(test_path);
+    cout << test_path << endl;
+    Grid grid = setupSimulation();
+    CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true);
+}
diff --git a/test/testSimulation.hpp.in b/test/testSimulation.hpp.in
new file mode 100644
index 0000000..2e1e19a
--- /dev/null
+++ b/test/testSimulation.hpp.in
@@ -0,0 +1,7 @@
+#ifndef TESTSIMULATION_H_
+#define TESTSIMULATION_H_
+
+// CSV file needed for validation
+const char *testSimulationCSVDir = "@testSimulationCSV@";
+
+#endif // TESTSIMULATION_H_