diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp
index c0936c4..3ac0df6 100644
--- a/examples/FTCS_1D_proto_example.cpp
+++ b/examples/FTCS_1D_proto_example.cpp
@@ -1,3 +1,4 @@
+#include "tug/Boundary.hpp"
 #include <tug/Simulation.hpp>
 
 int main(int argc, char *argv[]) {
@@ -18,7 +19,9 @@ int main(int argc, char *argv[]) {
     // ******************
 
     // create a boundary with constant values
-    Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
+    Boundary bc = Boundary(grid);
+    bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
+    bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
 
 
     // ************************
diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp
index 7171f62..deb778d 100644
--- a/examples/FTCS_2D_proto_example.cpp
+++ b/examples/FTCS_2D_proto_example.cpp
@@ -44,7 +44,12 @@ int main(int argc, char *argv[]) {
     // ******************
 
     // create a boundary with constant values
-    Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
+    Boundary bc = Boundary(grid);
+    bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
+    bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+    bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
+    bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
+
 
     // (optional) set boundary condition values for one side, e.g.:
     // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
@@ -68,7 +73,7 @@ int main(int argc, char *argv[]) {
     simulation.setTimestep(0.1); // timestep
 
     // (optional) set the number of iterations
-    simulation.setIterations(2);
+    simulation.setIterations(10000);
 
     // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
     simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
diff --git a/src/Boundary.cpp b/src/Boundary.cpp
index 9a1012a..87752ea 100644
--- a/src/Boundary.cpp
+++ b/src/Boundary.cpp
@@ -34,12 +34,15 @@ double BoundaryElement::getValue() {
 
 Boundary::Boundary(Grid grid) : grid(grid) {
     //probably to DEBUG assignment grid
-    
 
     if (grid.getDim() == 1) {
+        this->boundaries = vector<vector<BoundaryElement>>(2);
+
         this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement());
         this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement());
     } else if (grid.getDim() == 2) {
+        this->boundaries = vector<vector<BoundaryElement>>(4);
+
         this->boundaries[BC_SIDE_LEFT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
         this->boundaries[BC_SIDE_RIGHT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
         this->boundaries[BC_SIDE_TOP] = vector<BoundaryElement>(grid.getCol(), BoundaryElement());
diff --git a/src/FTCS.cpp b/src/FTCS.cpp
index 1b24e5a..f00891a 100644
--- a/src/FTCS.cpp
+++ b/src/FTCS.cpp
@@ -1,9 +1,11 @@
+#include "TugUtils.hpp"
 #include <cstddef>
 #include <tug/Boundary.hpp>
 #include <iostream>
 
 using namespace std;
 
+
 double calcAlphaIntercell(double alpha1, double alpha2, bool useHarmonic = false) {
     if (useHarmonic) {
         return 2 / ((1/alpha1) + (1/alpha2));
@@ -57,7 +59,7 @@ double calcHorizontalChangeLeftBoundaryConstant(Grid grid, Boundary bc, int row,
             + 2 * grid.getAlphaX()(row,col)
             ) 
             * grid.getConcentrations()(row,col) 
-        + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row);
+        + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
 
     return result;
 }
@@ -73,10 +75,21 @@ double calcHorizontalChangeLeftBoundaryClosed(Grid grid, int row, int col) {
 }
 
 
+double calcHorizontalChangeLeftBoundary(Grid grid, Boundary bc, int row, int col) {
+    if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) {
+        return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
+    } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) {
+        return calcHorizontalChangeLeftBoundaryClosed(grid, row, col);
+    } else {
+        throw_invalid_argument("Undefined Boundary Condition Type!");
+    }
+}
+
+
 double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
 
     double result = 
-        2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+        2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row)
         - (
             calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
             + 2 * grid.getAlphaX()(row,col)
@@ -99,6 +112,17 @@ double calcHorizontalChangeRightBoundaryClosed(Grid grid, int row, int col) {
 }
 
 
+double calcHorizontalChangeRightBoundary(Grid grid, Boundary bc, int row, int col) {
+    if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) {
+        return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
+    } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) {
+        return calcHorizontalChangeRightBoundaryClosed(grid, row, col);
+    } else {
+        throw_invalid_argument("Undefined Boundary Condition Type!");
+    }
+}
+
+
 double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
     
     double result = 
@@ -109,7 +133,7 @@ double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, in
             + 2 * grid.getAlphaY()(row, col)
             ) 
             * grid.getConcentrations()(row, col)
-        + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col);
+        + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col);
 
     return result;
 }
@@ -125,10 +149,21 @@ double calcVerticalChangeTopBoundaryClosed(Grid grid, int row, int col) {
 }
 
 
+double calcVerticalChangeTopBoundary(Grid grid, Boundary bc, int row, int col) {
+    if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
+        return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
+    } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
+        return calcVerticalChangeTopBoundaryClosed(grid, row, col);
+    } else {
+        throw_invalid_argument("Undefined Boundary Condition Type!");
+    }
+}
+
+
 double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
 
     double result = 
-        2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+        2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col)
         - (
             calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) 
             + 2 * grid.getAlphaY()(row, col)
@@ -151,6 +186,17 @@ double calcVerticalChangeBottomBoundaryClosed(Grid grid, int row, int col) {
 }
 
 
+double calcVerticalChangeBottomBoundary(Grid grid, Boundary bc, int row, int col) {
+    if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
+        return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
+    } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
+        return calcVerticalChangeBottomBoundaryClosed(grid, row, col);
+    } else {
+        throw_invalid_argument("Undefined Boundary Condition Type!");
+    }
+}
+
+
 MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
     int colMax = grid.getCol();
     double deltaCol = grid.getDeltaCol();
@@ -175,7 +221,7 @@ MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
     concentrations_t1(row, col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
                 * (
-                    calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
+                    calcHorizontalChangeLeftBoundary(grid, bc, row, col)
                 )
             ;
 
@@ -185,7 +231,7 @@ MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
                 * (
-                    calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
+                    calcHorizontalChangeRightBoundary(grid, bc, row, col)
                 )
             ;
 
@@ -240,7 +286,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
         concentrations_t1(row, col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
                 * (
-                    calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
+                    calcHorizontalChangeLeftBoundary(grid, bc, row, col)
                 )
             + timestep / (deltaRow*deltaRow)
                 * (
@@ -255,7 +301,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
         concentrations_t1(row,col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
                 * (
-                    calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
+                    calcHorizontalChangeRightBoundary(grid, bc, row, col)
                 )
             + timestep / (deltaRow*deltaRow)
                 * (
@@ -271,7 +317,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
             + timestep / (deltaRow*deltaRow) 
                 * (
-                    calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
+                    calcVerticalChangeTopBoundary(grid, bc, row, col)
                 )
             + timestep / (deltaCol*deltaCol) 
                 * (
@@ -286,7 +332,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
             + timestep / (deltaRow*deltaRow) 
                 * (
-                    calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
+                    calcVerticalChangeBottomBoundary(grid, bc, row, col)
                 )
             + timestep / (deltaCol*deltaCol) 
                 * (
@@ -301,11 +347,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
         + timestep/(deltaCol*deltaCol)
             * (
-                calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
+                calcHorizontalChangeLeftBoundary(grid, bc, row, col)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
+                calcVerticalChangeTopBoundary(grid, bc, row, col)
             )
         ;
 
@@ -315,11 +361,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
     concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
         + timestep/(deltaCol*deltaCol)
             * (
-                calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
+                calcHorizontalChangeRightBoundary(grid, bc, row, col)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
+                calcVerticalChangeTopBoundary(grid, bc, row, col)
             )
         ;
 
@@ -329,11 +375,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
         + timestep/(deltaCol*deltaCol)
             * (
-                calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
+                calcHorizontalChangeLeftBoundary(grid, bc, row, col)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
+                calcVerticalChangeBottomBoundary(grid, bc, row, col)
             )
         ;
 
@@ -343,11 +389,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
     concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
         + timestep/(deltaCol*deltaCol)
             * (
-                calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
+                calcHorizontalChangeRightBoundary(grid, bc, row, col)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
+                calcVerticalChangeBottomBoundary(grid, bc, row, col)
             )
         ;