diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp
index 9586454..c599448 100644
--- a/examples/FTCS_2D_proto_example.cpp
+++ b/examples/FTCS_2D_proto_example.cpp
@@ -1,16 +1,65 @@
+/**
+ * @file FTCS_2D_proto_example.cpp
+ * @author Hannes Signer, Philipp Ungrund
+ * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach
+ * and constant boundary condition
+ * 
+ */
+
 #include <tug/Simulation.hpp>
-#include <iostream>
 
 int main(int argc, char *argv[]) {
     
-    Grid grid = Grid(20,20);
+    // **************
+    // **** GRID ****
+    // **************
 
+    // create a grid with a 20 x 20 field
+    int row = 20;
+    int col = 20;
+    Grid grid = Grid(row,col);
+
+    // (optional) set the domain, e.g.:
+    // grid.setDomain(20, 20);
+
+    // (optional) set the concentrations, e.g.:
+    // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value
+    // grid.setConcentrations(concentrations);
+
+    // (optional) set alphas of the grid, e.g.:
+    // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value
+    // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value
+    // grid.setAlpha(alphax, alphay);
+
+
+    // ******************
+    // **** BOUNDARY ****
+    // ******************
+
+    // create a boundary with constant values
     Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
 
-    Simulation simulation = Simulation(grid, bc, FTCS_APPROACH);
-    simulation.setIterations(2);
+    // (optional) set boundary condition values for one side, e.g.:
+    // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
+    // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values
+
+
+    // ************************
+    // **** SIMULATION ENV ****
+    // ************************
+
+    // set up a simulation environment
+    Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
+
+    // (optional) set the timestep of the simulation
+    // simulation.setTimestep(0.01); // timestep
+
+    // (optional) set the number of iterations
+    simulation.setIterations(20);
     
+    // **** RUN SIMULATION ****
+
+    // run the simulation
     simulation.run();
-    
 
 }
\ No newline at end of file
diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp
index 3da1ee7..012e05d 100644
--- a/include/tug/Grid.hpp
+++ b/include/tug/Grid.hpp
@@ -20,6 +20,8 @@ class Grid {
         */
         Grid(int row, int col);
 
+        Grid(MatrixXd concentrations);
+
         /**
         * @brief Set the Concentrations object
         * 
diff --git a/src/FTCS.cpp b/src/FTCS.cpp
index 65cce82..6db41de 100644
--- a/src/FTCS.cpp
+++ b/src/FTCS.cpp
@@ -4,7 +4,7 @@
 
 using namespace std;
 
-auto calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) {
+double calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) {
     if (useHarmonic) {
         return 2 / ((1/alpha1) + (1/alpha2));
     } else {
@@ -23,14 +23,12 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
     MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 1);
 
     // inner cells
-    cout << "Concentration 5,5: " << grid.getConcentrations()(5,5) << endl;
-    cout << "Alpha Y 5,5: " << grid.getAlphaY()(5,5) << endl;
-    cout << "calc alpha Y 5,5; 5,6: " << calc_alpha_intercell(grid.getAlphaY()(5,5), grid.getAlphaY()(5,6)) << endl;
-    cout << "t1 Concentrations 5,5: " << concentrations_t1(5,5) << endl;
+    // (should have 7 calls to current concentration)
     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))
@@ -38,8 +36,9 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
                     * grid.getConcentrations()(row,col)
                     + 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))
@@ -47,94 +46,211 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
                     * grid.getConcentrations()(row,col)
                     + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
                     * grid.getConcentrations()(row,col-1)
-                );
+                    )
+                ;
         }
     }
 
     // boundary conditions
     // left without corners / looping over rows
+    // (should have 6 calls to current concentration)
     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))
-    //         + 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));
-    // }
+    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))
+            + 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));
+    }
 
     // right without corners / looping over columns
+    // (should have 6 calls to current concentration)
     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)
-    //             - (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));
-    // }
+    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)
+                - (calc_alpha_intercell(grid.getAlphaX()(row,col-1), 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)(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));
-    // }
+    // (should have 6 calls to current concentration)
+    int row = 0;
+    for (int col=1; col<colMax-1;col++){
+        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)(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)(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));
-    // }
+    // (should have 6 calls to current concentration)
+    row = rowMax-1;
+    for(int col=1; col<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)(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)) * grid.getConcentrations()(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)) * grid.getConcentrations()(row, col-1));
+    }
 
+    // corner top left
+    // (should have 5 calls to current concentration)
+    row = 0;
+    col = 0;
+    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)
+            )
+        + 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))
+                    + 2 * grid.getAlphaY()(row,col)
+                )
+                * grid.getConcentrations()(row,col)
+                + 2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col)
+            )
+        ;
+
+    // corner top right
+    // (should have 5 calls to current concentration)
+    row = 0;
+    col = colMax-1;
+    concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
+        + timestep/(deltaCol*deltaCol)
+            * (
+                2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                - (
+                    calc_alpha_intercell(grid.getAlphaX()(row,col-1), 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))
+                    + 2 * grid.getAlphaY()(row,col)
+                )
+                * grid.getConcentrations()(row,col)
+                + 2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col)
+            )
+        ;
+
+    // corner bottom left
+    // (should have 5 calls to current concentration)
+    row = rowMax-1;
+    col = 0;
+    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)
+            )
+        + timestep/(deltaRow*deltaRow)
+            * (
+                2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(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))
+                * grid.getConcentrations()(row-1,col)
+            )
+        ;
+
+    // corner bottom right
+    // (should have 5 calls to current concentration)
+    row = rowMax-1;
+    col = colMax-1;
+    concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
+        + timestep/(deltaCol*deltaCol)
+            * (
+                2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                - (
+                    calc_alpha_intercell(grid.getAlphaX()(row,col-1), 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)
+            * (
+                2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(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))
+                * grid.getConcentrations()(row-1,col)
+            )
+        ;
 
-    concentrations_t1(0,0) = 0;
-    concentrations_t1(rowMax-1,0) = 0;
-    concentrations_t1(0,colMax-1) = 0;
-    concentrations_t1(rowMax-1,colMax-1) = 0;
 
     return concentrations_t1;
 }
 
-void FTCS_closed(Grid grid, Boundary bc, double timestep) {
-    return;
+// TODO
+MatrixXd FTCS_closed(Grid grid, Boundary bc, double timestep) {
+    return MatrixXd();
 }
 
 MatrixXd FTCS(Grid grid, Boundary bc, double timestep) {
-    if (bc.getBoundaryConditionType() == BC_TYPE_CONSTANT) {
-        return FTCS_constant(grid, bc, timestep);
-    } else if (bc.getBoundaryConditionType() == BC_TYPE_CLOSED) {
-        FTCS_closed(grid, bc, timestep);
+    switch (bc.getBoundaryConditionType()) {
+        case BC_TYPE_CONSTANT:
+            return FTCS_constant(grid, bc, timestep);
+        case BC_TYPE_CLOSED:
+            return FTCS_closed(grid, bc, timestep);
+        default:
+            // TODO handle
+            return MatrixXd();
     }
 }
diff --git a/src/Grid.cpp b/src/Grid.cpp
index f00f720..28b42f7 100644
--- a/src/Grid.cpp
+++ b/src/Grid.cpp
@@ -3,18 +3,24 @@
 
 Grid::Grid(int col) {
     this->col = col;
+    this->domain_col = col;
+    this->delta_col = double(this->domain_col)/double(this->col);
 
     this->dim = 1;
-    this->concentrations = MatrixXd::Constant(1, col, 1);
+    this->concentrations = MatrixXd::Constant(1, col, 20);
     this->alpha_x = MatrixXd::Constant(1, col, 1);
 }
 
 Grid::Grid(int row, int col) {
     this->row = row;
     this->col = col;
+    this->domain_row = row;
+    this->domain_col = col;
+    this->delta_row = double(this->domain_row)/double(this->row);
+    this->delta_col = double(this->domain_col)/double(this->col);
 
     this->dim = 2;
-    this->concentrations = MatrixXd::Constant(row, col, 1);
+    this->concentrations = MatrixXd::Constant(row, col, 20);
     this->alpha_x = MatrixXd::Constant(row, col, 1);
     this->alpha_y = MatrixXd::Constant(row, col, 1);
     
@@ -66,8 +72,8 @@ void Grid::setDomain(int domain_row, int domain_col) {
     this->domain_row = domain_row;
     this->domain_col = domain_col;
 
-    this->delta_row = double(this->domain_row)/this->row;
-    this->delta_col = double(this->domain_col)/this->col;
+    this->delta_row = double(this->domain_row)/double(this->row);
+    this->delta_col = double(this->domain_col)/double(this->col);
 }
 
 double Grid::getDeltaCol() {