diff --git a/julia/tug/Core/FTCS.jl b/julia/tug/Core/FTCS.jl
index 87d9fd5..4afd953 100644
--- a/julia/tug/Core/FTCS.jl
+++ b/julia/tug/Core/FTCS.jl
@@ -7,152 +7,88 @@ include("../Boundary.jl")
 include("../Grid.jl")
 include("Utils.jl")
 
-function calcHorizontalChange(grid::Grid{T}, row::Int, col::Int) where {T}
-    return calcAlphaIntercell(getAlphaX(grid)[row, col+1],
-        getAlphaX(grid)[row, col]) *
-           getConcentrations(grid)[row, col+1] -
-           (calcAlphaIntercell(getAlphaX(grid)[row, col+1],
-        getAlphaX(grid)[row, col]) +
-            calcAlphaIntercell(getAlphaX(grid)[row, col-1],
-        getAlphaX(grid)[row, col])) *
-           getConcentrations(grid)[row, col] +
-           calcAlphaIntercell(getAlphaX(grid)[row, col-1],
-        getAlphaX(grid)[row, col]) *
-           getConcentrations(grid)[row, col-1]
+function calcHorizontalChange(alphaX_next::T, alphaX_prev::T, alphaX_current::T,
+    concentration_next::T, concentration_prev::T, concentration_current::T) where {T}
+
+    intercellAlpha_next = calcAlphaIntercell(alphaX_next, alphaX_current)
+    intercellAlpha_prev = calcAlphaIntercell(alphaX_prev, alphaX_current)
+
+    return intercellAlpha_next * concentration_next -
+           (intercellAlpha_next + intercellAlpha_prev) * concentration_current +
+           intercellAlpha_prev * concentration_prev
 end
 
-function calcHorizontalChangeLeftBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    if getBoundaryElementType(bc, LEFT, row) == CONSTANT
-        return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
-    elseif getBoundaryElementType(bc, LEFT, row) == CLOSED
-        return calcHorizontalChangeLeftBoundaryClosed(grid, row, col)
+function calcHorizontalChangeLeftBoundary(boundaryType::TYPE, alphaX_next::T, alphaX_current::T,
+    concentration_next::T, concentration_current::T, boundaryValue::T) where {T}
+    intercellAlpha = calcAlphaIntercell(alphaX_next, alphaX_current)
+
+    if boundaryType == CONSTANT
+        return intercellAlpha * concentration_next -
+               (intercellAlpha + 2 * alphaX_current) * concentration_current +
+               2 * alphaX_current * boundaryValue
+    elseif boundaryType == CLOSED
+        return intercellAlpha * (concentration_next - concentration_current)
     else
         error("Undefined Boundary Condition Type!")
     end
 end
 
-function calcHorizontalChangeLeftBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    return calcAlphaIntercell(getAlphaX(grid)[row, col+1],
-        getAlphaX(grid)[row, col]) *
-           getConcentrations(grid)[row, col+1] -
-           (calcAlphaIntercell(getAlphaX(grid)[row, col+1],
-        getAlphaX(grid)[row, col]) +
-            2 * getAlphaX(grid)[row, col]) *
-           getConcentrations(grid)[row, col] +
-           2 * getAlphaX(grid)[row, col] *
-           getBoundaryElementValue(bc, LEFT, row)
-end
+function calcHorizontalChangeRightBoundary(boundaryType::TYPE, alphaX_prev::T, alphaX_current::T,
+    concentration_prev::T, concentration_current::T, boundaryValue::T) where {T}
+    intercellAlpha = calcAlphaIntercell(alphaX_prev, alphaX_current)
 
-function calcHorizontalChangeLeftBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T}
-    return calcAlphaIntercell(getAlphaX(grid)[row, col+1],
-        getAlphaX(grid)[row, col]) *
-           (getConcentrations(grid)[row, col+1] -
-            getConcentrations(grid)[row, col])
-end
-
-
-function calcHorizontalChangeRightBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    if getBoundaryElementType(bc, RIGHT, row) == CONSTANT
-        return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
-    elseif getBoundaryElementType(bc, RIGHT, row) == CLOSED
-        return calcHorizontalChangeRightBoundaryClosed(grid, row, col)
+    if boundaryType == CONSTANT
+        return 2 * alphaX_current * boundaryValue -
+               (intercellAlpha + 2 * alphaX_current) * concentration_current +
+               intercellAlpha * concentration_prev
+    elseif boundaryType == CLOSED
+        return -intercellAlpha * (concentration_current - concentration_prev)
     else
         error("Undefined Boundary Condition Type!")
     end
 end
 
-function calcHorizontalChangeRightBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    return 2 * getAlphaX(grid)[row, col] *
-           getBoundaryElementValue(bc, RIGHT, row) -
-           (calcAlphaIntercell(getAlphaX(grid)[row, col-1],
-        getAlphaX(grid)[row, col]) +
-            2 * getAlphaX(grid)[row, col]) *
-           getConcentrations(grid)[row, col] +
-           calcAlphaIntercell(getAlphaX(grid)[row, col-1],
-        getAlphaX(grid)[row, col]) *
-           getConcentrations(grid)[row, col-1]
+function calcVerticalChange(alphaY_above::T, alphaY_below::T, alphaY_current::T,
+    concentration_above::T, concentration_below::T, concentration_current::T) where {T}
+    intercellAlpha_above = calcAlphaIntercell(alphaY_above, alphaY_current)
+    intercellAlpha_below = calcAlphaIntercell(alphaY_below, alphaY_current)
+
+    return intercellAlpha_above * concentration_above -
+           (intercellAlpha_above + intercellAlpha_below) * concentration_current +
+           intercellAlpha_below * concentration_below
 end
 
-function calcHorizontalChangeRightBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T}
-    return -(calcAlphaIntercell(getAlphaX(grid)[row, col-1],
-        getAlphaX(grid)[row, col]) *
-             (getConcentrations(grid)[row, col] -
-              getConcentrations(grid)[row, col-1]))
-end
+function calcVerticalChangeTopBoundary(boundaryType::TYPE, alphaY_above::T, alphaY_current::T,
+    concentration_above::T, concentration_current::T, boundaryValue::T) where {T}
+    intercellAlpha = calcAlphaIntercell(alphaY_above, alphaY_current)
 
-
-function calcVerticalChange(grid::Grid{T}, row::Int, col::Int) where {T}
-    return calcAlphaIntercell(getAlphaY(grid)[row+1, col],
-        getAlphaY(grid)[row, col]) *
-           getConcentrations(grid)[row+1, col] -
-           (calcAlphaIntercell(getAlphaY(grid)[row+1, col],
-        getAlphaY(grid)[row, col]) +
-            calcAlphaIntercell(getAlphaY(grid)[row-1, col],
-        getAlphaY(grid)[row, col])) *
-           getConcentrations(grid)[row, col] +
-           calcAlphaIntercell(getAlphaY(grid)[row-1, col],
-        getAlphaY(grid)[row, col]) *
-           getConcentrations(grid)[row-1, col]
-end
-
-function calcVerticalChangeTopBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    if getBoundaryElementType(bc, TOP, col) == CONSTANT
-        return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
-    elseif getBoundaryElementType(bc, TOP, col) == CLOSED
-        return calcVerticalChangeTopBoundaryClosed(grid, row, col)
+    if boundaryType == CONSTANT
+        return intercellAlpha * concentration_above -
+               (intercellAlpha + 2 * alphaY_current) * concentration_current +
+               2 * alphaY_current * boundaryValue
+    elseif boundaryType == CLOSED
+        return intercellAlpha * (concentration_above - concentration_current)
     else
         error("Undefined Boundary Condition Type!")
     end
 end
 
-function calcVerticalChangeTopBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    return calcAlphaIntercell(getAlphaY(grid)[row+1, col],
-        getAlphaY(grid)[row, col]) *
-           getConcentrations(grid)[row+1, col] -
-           (calcAlphaIntercell(getAlphaY(grid)[row+1, col],
-        getAlphaY(grid)[row, col]) +
-            2 * getAlphaY(grid)[row, col]) *
-           getConcentrations(grid)[row, col] +
-           2 * getAlphaY(grid)[row, col] *
-           getBoundaryElementValue(bc, TOP, col)
-end
 
-function calcVerticalChangeTopBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T}
-    return calcAlphaIntercell(getAlphaY(grid)[row+1, col],
-        getAlphaY(grid)[row, col]) *
-           (getConcentrations(grid)[row+1, col] -
-            getConcentrations(grid)[row, col])
-end
+function calcVerticalChangeBottomBoundary(boundaryType::TYPE, alphaY_current::T, alphaY_below::T,
+    concentration_current::T, concentration_below::T, boundaryValue::T) where {T}
+    intercellAlpha = calcAlphaIntercell(alphaY_current, alphaY_below)
 
-function calcVerticalChangeBottomBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    if getBoundaryElementType(bc, BOTTOM, col) == CONSTANT
-        return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
-    elseif getBoundaryElementType(bc, BOTTOM, col) == CLOSED
-        return calcVerticalChangeBottomBoundaryClosed(grid, row, col)
+    if boundaryType == CONSTANT
+        return 2 * alphaY_current * boundaryValue -
+               (intercellAlpha + 2 * alphaY_current) * concentration_current +
+               intercellAlpha * concentration_below
+    elseif boundaryType == CLOSED
+        return -intercellAlpha * (concentration_current - concentration_below)
     else
         error("Undefined Boundary Condition Type!")
     end
 end
 
-function calcVerticalChangeBottomBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T}
-    return 2 * getAlphaY(grid)[row, col] *
-           getBoundaryElementValue(bc, BOTTOM, col) -
-           (calcAlphaIntercell(getAlphaY(grid)[row, col],
-        getAlphaY(grid)[row-1, col]) +
-            2 * getAlphaY(grid)[row, col]) *
-           getConcentrations(grid)[row, col] +
-           calcAlphaIntercell(getAlphaY(grid)[row, col],
-        getAlphaY(grid)[row-1, col]) *
-           getConcentrations(grid)[row-1, col]
-end
-
-function calcVerticalChangeBottomBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T}
-    return -(calcAlphaIntercell(getAlphaY(grid)[row, col],
-        getAlphaY(grid)[row-1, col]) *
-             (getConcentrations(grid)[row, col] -
-              getConcentrations(grid)[row-1, col]))
-end
-
 
 
 function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
@@ -160,78 +96,119 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
     sx = timestep / (getDeltaCol(grid)^2)
 
     concentrations = getConcentrations(grid)
+    alphaX = getAlphaX(grid)
     concentrations_t1 = copy(concentrations)
 
-    # only one row in 1D case -> row constant at index 1
     row = 1
 
     # inner cells
-    # independent of boundary condition type
     for col = 2:colMax-1
-        concentrations_t1[row, col] += sx * calcHorizontalChange(grid, row, col)
+        concentrations_t1[row, col] += sx * calcHorizontalChange(alphaX[row, col+1], alphaX[row, col-1], alphaX[row, col],
+            concentrations[row, col+1], concentrations[row, col-1], concentrations[row, col])
     end
 
-    # left boundary; hold column constant at index 1
-    concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, row, 1)
+    # left boundary
+    leftBoundaryType = getBoundaryElementType(bc, LEFT, row)
+    leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row)
+    concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(leftBoundaryType, alphaX[row, 2], alphaX[row, 1],
+        concentrations[row, 2], concentrations[row, 1], leftBoundaryValue)
 
-    # right boundary; hold column constant at max index
-    concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, row, colMax)
+    # right boundary
+    rightBoundaryType = getBoundaryElementType(bc, RIGHT, row)
+    rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row)
+    concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(rightBoundaryType, alphaX[row, colMax-1], alphaX[row, colMax],
+        concentrations[row, colMax-1], concentrations[row, colMax], rightBoundaryValue)
 
-    # overwrite obsolete concentrations
     setConcentrations!(grid, concentrations_t1)
 end
 
+
 function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
     rowMax = getRows(grid)
     colMax = getCols(grid)
     sx = timestep / (getDeltaCol(grid)^2)
     sy = timestep / (getDeltaRow(grid)^2)
 
+    alphaX = getAlphaX(grid)
+    alphaY = getAlphaY(grid)
     concentrations = getConcentrations(grid)
     concentrations_t1 = copy(concentrations)
 
+    # Currently boundary sides can only be set for the whole side, not for each cell individually
+    leftBoundaryType = getBoundaryElementType(bc, LEFT, 1)
+    rightBoundaryType = getBoundaryElementType(bc, RIGHT, 1)
+    topBoundaryType = getBoundaryElementType(bc, TOP, 1)
+    bottomBoundaryType = getBoundaryElementType(bc, BOTTOM, 1)
+
     Threads.@threads for row = 2:rowMax-1
-        # inner cells
         for col = 2:colMax-1
-            concentrations_t1[row, col] += sy * calcVerticalChange(grid, row, col) +
-                                           sx * calcHorizontalChange(grid, row, col)
+            concentrations_t1[row, col] += sy * calcVerticalChange(alphaY[row+1, col], alphaY[row-1, col], alphaY[row, col],
+                concentrations[row+1, col], concentrations[row-1, col], concentrations[row, col]) +
+                                           sx * calcHorizontalChange(alphaX[row, col+1], alphaX[row, col-1], alphaX[row, col],
+                concentrations[row, col+1], concentrations[row, col-1], concentrations[row, col])
         end
 
-        # left boundary without corners
-        concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, row, 1) +
-                                     sy * calcVerticalChange(grid, row, 1)
+        # Boundary conditions for each row
+        # Left boundary without corners
+        leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row)
+        concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(leftBoundaryType, alphaX[row, 2], alphaX[row, 1],
+            concentrations[row, 2], concentrations[row, 1], leftBoundaryValue) +
+                                     sy * calcVerticalChange(alphaY[row+1, 1], alphaY[row-1, 1], alphaY[row, 1],
+            concentrations[row+1, 1], concentrations[row-1, 1], concentrations[row, 1])
 
-        # right boundary without corners
-        concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, row, colMax) +
-                                          sy * calcVerticalChange(grid, row, colMax)
+        # Right boundary without corners
+        rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row)
+        concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(rightBoundaryType, alphaX[row, colMax-1], alphaX[row, colMax],
+            concentrations[row, colMax-1], concentrations[row, colMax], rightBoundaryValue) +
+                                          sy * calcVerticalChange(alphaY[row+1, colMax], alphaY[row-1, colMax], alphaY[row, colMax],
+            concentrations[row+1, colMax], concentrations[row-1, colMax], concentrations[row, colMax])
     end
 
-    # top / bottom without corners / looping over columns
+    # Handle top/bottom boundaries
     Threads.@threads for col = 2:colMax-1
-        # top
-        concentrations_t1[1, col] += sy * calcVerticalChangeTopBoundary(grid, bc, 1, col) +
-                                     sx * calcHorizontalChange(grid, 1, col)
+        # Top boundary
+        topBoundaryValue = getBoundaryElementValue(bc, TOP, col)
+        concentrations_t1[1, col] += sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, col], alphaY[1, col],
+            concentrations[2, col], concentrations[1, col], topBoundaryValue) +
+                                     sx * calcHorizontalChange(alphaX[1, col+1], alphaX[1, col-1], alphaX[1, col],
+            concentrations[1, col+1], concentrations[1, col-1], concentrations[1, col])
 
-        # bottom
-        concentrations_t1[rowMax, col] += sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, col) +
-                                          sx * calcHorizontalChange(grid, rowMax, col)
+        # Bottom boundary
+        bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, col)
+        concentrations_t1[rowMax, col] += sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, col], alphaY[rowMax-1, col],
+            concentrations[rowMax, col], concentrations[rowMax-1, col], bottomBoundaryValue) +
+                                          sx * calcHorizontalChange(alphaX[rowMax, col+1], alphaX[rowMax, col-1], alphaX[rowMax, col],
+            concentrations[rowMax, col+1], concentrations[rowMax, col-1], concentrations[rowMax, col])
     end
 
-    # corner top left
-    concentrations_t1[1, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, 1, 1) +
-                               sy * calcVerticalChangeTopBoundary(grid, bc, 1, 1)
+    # Handle corners
+    # Top left corner
+    topBoundaryValue = getBoundaryElementValue(bc, TOP, 1)
+    concentrations_t1[1, 1] += sx * calcHorizontalChange(alphaX[1, 2], alphaX[1, 1], alphaX[1, 1],
+        concentrations[1, 2], concentrations[1, 1], concentrations[1, 1]) +
+                               sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, 1], alphaY[1, 1],
+        concentrations[2, 1], concentrations[1, 1], topBoundaryValue)
 
-    # corner top right
-    concentrations_t1[1, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, 1, colMax) +
-                                    sy * calcVerticalChangeTopBoundary(grid, bc, 1, colMax)
+    # Top right corner
+    topBoundaryValue = getBoundaryElementValue(bc, TOP, colMax)
+    concentrations_t1[1, colMax] += sx * calcHorizontalChange(alphaX[1, colMax-1], alphaX[1, colMax], alphaX[1, colMax],
+        concentrations[1, colMax-1], concentrations[1, colMax], concentrations[1, colMax]) +
+                                    sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, colMax], alphaY[1, colMax],
+        concentrations[2, colMax], concentrations[1, colMax], topBoundaryValue)
 
-    # corner bottom left
-    concentrations_t1[rowMax, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, rowMax, 1) +
-                                    sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, 1)
+    # Bottom left corner
+    bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, 1)
+    concentrations_t1[rowMax, 1] += sx * calcHorizontalChange(alphaX[rowMax, 2], alphaX[rowMax, 1], alphaX[rowMax, 1],
+        concentrations[rowMax, 2], concentrations[rowMax, 1], concentrations[rowMax, 1]) +
+                                    sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, 1], alphaY[rowMax-1, 1],
+        concentrations[rowMax, 1], concentrations[rowMax-1, 1], bottomBoundaryValue)
 
-    # corner bottom right
-    concentrations_t1[rowMax, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, rowMax, colMax) +
-                                         sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, colMax)
+    # Bottom right corner
+    bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, colMax)
+    concentrations_t1[rowMax, colMax] += sx * calcHorizontalChange(alphaX[rowMax, colMax-1], alphaX[rowMax, colMax], alphaX[rowMax, colMax],
+        concentrations[rowMax, colMax-1], concentrations[rowMax, colMax], concentrations[rowMax, colMax]) +
+                                         sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, colMax], alphaY[rowMax-1, colMax],
+        concentrations[rowMax, colMax], concentrations[rowMax-1, colMax], bottomBoundaryValue)
 
     setConcentrations!(grid, concentrations_t1)
 end