diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 9723c76..d46ed8f 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,7 +1,9 @@
 add_executable(first_example first_example.cpp)
 add_executable(second_example second_example.cpp)
 add_executable(boundary_example1D boundary_example1D.cpp)
+add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
 
 target_link_libraries(first_example tug)
 target_link_libraries(second_example tug)
-target_link_libraries(boundary_example1D tug)
\ No newline at end of file
+target_link_libraries(boundary_example1D tug)
+target_link_libraries(FTCS_2D_proto_example tug)
\ No newline at end of file
diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp
new file mode 100644
index 0000000..2819334
--- /dev/null
+++ b/examples/FTCS_2D_proto_example.cpp
@@ -0,0 +1,11 @@
+#include <tug/Simulation.hpp>
+
+int main(int argc, char *argv[]) {
+    Grid grid = Grid(20,20);
+
+    Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
+
+    Simulation simulation = Simulation(grid, bc, FTCS_APPROACH);
+
+    simulation.run();
+}
\ No newline at end of file
diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp
index 383a6b9..c49d456 100644
--- a/include/tug/Boundary.hpp
+++ b/include/tug/Boundary.hpp
@@ -50,7 +50,7 @@ class Boundary {
          * @param side 
          * @return auto 
          */
-        auto getBoundaryConditionValue(BC_SIDE side);
+        VectorXd getBoundaryConditionValue(BC_SIDE side);
 
 
     private:
diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp
index cbb8c64..8d0d809 100644
--- a/include/tug/Simulation.hpp
+++ b/include/tug/Simulation.hpp
@@ -62,7 +62,7 @@ class Simulation {
          * 
          * @return auto 
          */
-        auto run();
+        void run();
 
     private:
 
diff --git a/src/Boundary.cpp b/src/Boundary.cpp
index 2c25848..03e5216 100644
--- a/src/Boundary.cpp
+++ b/src/Boundary.cpp
@@ -51,7 +51,7 @@ void Boundary::setBoundaryConditionValue(BC_SIDE side, VectorXd &values) {
     }
 }
 
-auto Boundary::getBoundaryConditionValue(BC_SIDE side) {
+VectorXd Boundary::getBoundaryConditionValue(BC_SIDE side) {
     switch (side) {
         case BC_SIDE_LEFT:
             return this->left;
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 35932a1..3790689 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,4 +1,4 @@
-add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp)
+add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp)
 
 target_link_libraries(tug Eigen3::Eigen)
 
diff --git a/src/FTCS.cpp b/src/FTCS.cpp
index f7aca0d..c32ca97 100644
--- a/src/FTCS.cpp
+++ b/src/FTCS.cpp
@@ -23,7 +23,7 @@ auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) {
     for (int row = 1; row < rowMax-1; row++) {
         for (int col = 1; col < colMax-1; col++) {
             concentrations_t1(row, col) = grid.getConcentrations()(row, col) 
-                + timestep / deltaRow*deltaRow * (
+                + timestep / (deltaRow*deltaRow) * (
                     calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
                     * grid.getConcentrations()(row+1,col)
                     - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
@@ -32,7 +32,7 @@ auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) {
                     + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
                     * grid.getConcentrations()(row-1,col)
                 )
-                - timestep / deltaCol*deltaCol * (
+                - timestep / (deltaCol*deltaCol) * (
                     calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
                     * grid.getConcentrations()(row,col+1)
                     - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
@@ -46,17 +46,87 @@ auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) {
 
     // boundary conditions
     // left without corners / looping over rows
+    int col = 0;
+    for (int row = 1; row < rowMax-1; row++) {
+        concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+            + timestep / (deltaCol*deltaCol) 
+                * (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                * grid.getConcentrations()(row,col+1)
+                - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                + 2 * grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col) 
+                + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row, 1))
+            + timestep / (deltaRow*deltaRow)
+                * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                * grid.getConcentrations()(row+1,col)
+                - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
+                * grid.getConcentrations()(row,col)
+                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getConcentrations()(row,col))
+                * grid.getConcentrations()(row-1,col));
+    }
 
-    return 0;
+    // right without corners / looping over columns
+    col = colMax-1;
+    for (int row = 1; row < rowMax-1; row++) {
+        concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+            + timestep / (deltaCol*deltaCol) 
+                * (2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row, 1)
+                - (calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+                + 2 * grid.getAlphaX()(row,col)) + 2 * grid.getAlphaX()(row,col) 
+                * grid.getConcentrations()(row,col)
+                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+                * grid.getConcentrations()(row,col-1))
+            + timestep / (deltaRow*deltaRow)
+                * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+                * grid.getConcentrations()(row+1,col)
+                - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
+                * grid.getConcentrations()(row,col)
+                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
+                * grid.getConcentrations()(row-1,col));
+    }
+
+
+    // top without corners / looping over cols
+    for(int col=1; col<colMax-1;col++){
+        int row = 0;
+        concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+            + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) * grid.getConcentrations()(1,col)
+                - (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) + 2 * grid.getAlphaY()(0, col)) * grid.getConcentrations()(0, col)
+                + 2 * grid.getAlphaY()(0, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(1, col))
+            + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col+1)
+                - (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col))) * grid.getConcentrations()(0, col)
+                + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col-1));
+    }
+
+    // bottom without corners / looping over cols
+    int row = rowMax-1;
+    for(int col=1; row<colMax-1;col++){
+        concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+            + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(1, col)
+                - (calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + 2 * grid.getAlphaY()(row, col)) * grid.getConcentrations()(row, col)
+                + calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)))
+            + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) * grid.getConcentrations()(row, col+1)
+                - (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))) * grid.getConcentrations()(row, col)
+                + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col-1)) * grid.getConcentrations()(row, col-1));
+    }
+
+
+    concentrations_t1(0,0) = 0;
+    concentrations_t1(rowMax-1,0) = 0;
+    concentrations_t1(0,colMax-1) = 0;
+    concentrations_t1(rowMax-1,colMax-1) = 0;
+
+    grid.setConcentrations(concentrations_t1);
 }
 
 auto FTCS_closed(Grid &grid, Boundary &bc, double timestep) {
     return 0;
 }
 
-auto FTCS(Grid &grid, Boundary &bc, double timestep) {
+void FTCS(Grid &grid, Boundary &bc, double timestep) {
     if (bc.getBoundaryConditionType() == BC_TYPE_CONSTANT) {
-        int test = FTCS_constant(grid, bc, timestep);
+        FTCS_constant(grid, bc, timestep);
     } else if (bc.getBoundaryConditionType() == BC_TYPE_CLOSED) {
         FTCS_closed(grid, bc, timestep);
     }
diff --git a/src/Simulation.cpp b/src/Simulation.cpp
index 0e44d32..92a042a 100644
--- a/src/Simulation.cpp
+++ b/src/Simulation.cpp
@@ -46,7 +46,7 @@ auto Simulation::getIterations() {
 }
 
 
-auto Simulation::run() {
+void Simulation::run() {
     if (approach == FTCS_APPROACH) {
         for (int i = 0; i < iterations; i++) {
             FTCS(grid, bc, timestep);