diff --git a/julia/TUG/src/AbstractSimulation.jl b/julia/TUG/src/AbstractSimulation.jl
index 69ebb0c..05b8e6b 100644
--- a/julia/TUG/src/AbstractSimulation.jl
+++ b/julia/TUG/src/AbstractSimulation.jl
@@ -1,7 +1,33 @@
 using Printf
 
+"""
+    @enum APPROACH BTCS FTCS
+
+Enum representing the numerical approach for the simulation.
+- `BTCS`: Backward Time, Central Space. A stable, implicit method.
+- `FTCS`: Forward Time, Central Space. A conditionally stable, explicit method.
+"""
 @enum APPROACH BTCS FTCS
+
+"""
+    @enum CONSOLE_OUTPUT CONSOLE_OUTPUT_OFF CONSOLE_OUTPUT_ON CONSOLE_OUTPUT_VERBOSE
+
+Enum representing the level of console output for the simulation.
+- `CONSOLE_OUTPUT_OFF`: No output to the console.
+- `CONSOLE_OUTPUT_ON`: Standard output level.
+- `CONSOLE_OUTPUT_VERBOSE`: Detailed output, including intermediate steps.
+"""
 @enum CONSOLE_OUTPUT CONSOLE_OUTPUT_OFF CONSOLE_OUTPUT_ON CONSOLE_OUTPUT_VERBOSE
+
+"""
+    @enum CSV_OUTPUT CSV_OUTPUT_OFF CSV_OUTPUT_ON CSV_OUTPUT_VERBOSE CSV_OUTPUT_XTREME
+
+Enum representing the level of CSV file output for the simulation.
+- `CSV_OUTPUT_OFF`: No output to CSV.
+- `CSV_OUTPUT_ON`: Standard output level.
+- `CSV_OUTPUT_VERBOSE`: Detailed output.
+- `CSV_OUTPUT_XTREME`: Extremely detailed output, possibly including debugging information.
+"""
 @enum CSV_OUTPUT CSV_OUTPUT_OFF CSV_OUTPUT_ON CSV_OUTPUT_VERBOSE CSV_OUTPUT_XTREME
 
 abstract type AbstractSimulation{T} end
@@ -38,7 +64,7 @@ function adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iteratio
     return timestep, iterations
 end
 
-function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing)::IOStream where {T}
+function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing) where {T}
     grid = grid === nothing ? simulation.grid : grid
     approachString = string(simulation.approach)
     rows = getRows(grid)
@@ -68,8 +94,7 @@ function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},No
         end
     end
 
-    file = open(filename, "a")
-    return file
+    open(filename, "a")
 end
 
 function writeBoundarySideValuesToCSV(file::IOStream, bc::Boundary{T}, side) where {T}
diff --git a/julia/TUG/src/Boundary.jl b/julia/TUG/src/Boundary.jl
index 3aa3e31..716784e 100644
--- a/julia/TUG/src/Boundary.jl
+++ b/julia/TUG/src/Boundary.jl
@@ -3,18 +3,49 @@
 # condition at the edges of the diffusion grid. 
 # Translated from C++'s Boundary.hpp.
 
+"""
+    @enum TYPE CLOSED CONSTANT
+
+Enum representing the type of boundary.
+- `CLOSED`: Closed boundary condition.
+- `CONSTANT`: Constant boundary condition.
+"""
 @enum TYPE CLOSED CONSTANT
+
+"""
+    @enum SIDE LEFT RIGHT TOP BOTTOM
+
+Enum representing the side of the boundary.
+- `LEFT`: Left side.
+- `RIGHT`: Right side.
+- `TOP`: Top side.
+- `BOTTOM`: Bottom side.
+"""
 @enum SIDE LEFT = 1 RIGHT = 2 TOP = 3 BOTTOM = 4
 
+"""
+    struct BoundaryElement{T}
+
+Represents an element of a boundary condition in a diffusion grid. 
+It holds the type of boundary condition and its value if applicable.
+
+# Fields
+- `type::TYPE`: The type of boundary condition (CLOSED, CONSTANT).
+- `value::T`: The value of the boundary condition, used if type is CONSTANT.
+
+# Constructors
+- `BoundaryElement{T}()` creates a CLOSED boundary element.
+- `BoundaryElement{T}(value::T)` creates a CONSTANT boundary element with specified value.
+"""
 struct BoundaryElement{T}
     type::TYPE
     value::T
 
-    # Generic constructor for closed case
-    BoundaryElement{T}() where {T} = new{T}(CLOSED, convert(T, -1))
+    function BoundaryElement{T}()::BoundaryElement{T} where {T}
+        new{T}(CLOSED, convert(T, -1))
+    end
 
-    # Constructor for constant case
-    function BoundaryElement{T}(value::T) where {T}
+    function BoundaryElement{T}(value::T)::BoundaryElement{T} where {T}
         if value < 0
             throw(ArgumentError("No negative concentration allowed."))
         end
@@ -22,26 +53,72 @@ struct BoundaryElement{T}
     end
 end
 
+"""
+    getType(be::BoundaryElement{T})::TYPE where {T}
+
+Retrieves the type of the boundary element.
+
+# Arguments
+- `be::BoundaryElement{T}`: The boundary element.
+
+# Returns
+The type of the boundary element (`TYPE`).
+"""
 function getType(be::BoundaryElement{T})::TYPE where {T}
     be.type
 end
 
+"""
+    getValue(be::BoundaryElement{T})::T where {T}
+
+Retrieves the value of the boundary element. Applicable if the boundary type is CONSTANT.
+
+# Arguments
+- `be::BoundaryElement{T}`: The boundary element.
+
+# Returns
+The value of the boundary element.
+"""
 function getValue(be::BoundaryElement{T})::T where {T}
     be.value
 end
 
+"""
+    getValue(be::Vector{BoundaryElement{T}})::Vector{T} where {T}
+
+Retrieves the values of the boundary elements. Applicable if the boundary type is CONSTANT.
+
+# Arguments
+- `be::Vector{BoundaryElement{T}}`: The boundary elements.
+
+# Returns
+The values of the boundary elements.
+"""
 function getValue(be::Vector{BoundaryElement{T}})::Vector{T} where {T}
     [b.value for b in be]
 end
 
+"""
+    struct Boundary{T}
+
+Represents the boundary conditions of a diffusion grid. 
+It includes the dimension of the grid and boundary elements for each side.
+
+# Fields
+- `dim::UInt8`: Dimension of the grid (1 or 2).
+- `cols::UInt32`, `rows::UInt32`: Number of columns and rows in the grid.
+- `boundaries::Vector{Vector{BoundaryElement{T}}}`: A vector of boundary elements for each side.
+
+# Constructor
+- `Boundary(grid::Grid{T})` constructs the boundary based on the given grid's dimensions.
+"""
 struct Boundary{T}
     dim::UInt8
     cols::UInt32
     rows::UInt32
     boundaries::Vector{Vector{BoundaryElement{T}}}
 
-    # Constructor
-    function Boundary(grid::Grid{T}) where {T}
+    function Boundary(grid::Grid{T})::Boundary{T} where {T}
         dim = grid.dim
         cols = grid.cols
         rows = grid.rows
@@ -63,6 +140,19 @@ struct Boundary{T}
     end
 end
 
+"""
+    getBoundaryElementType(boundary::Boundary{T}, side::SIDE, index::Int)::TYPE where {T}
+
+Retrieves the type of the boundary element at the specified side and index.
+
+# Arguments
+- `boundary::Boundary{T}`: The boundary.
+- `side::SIDE`: The side of the boundary.
+- `index::Int`: The index of the boundary element.
+
+# Returns
+The type of the boundary element (`TYPE`).
+"""
 function getBoundaryElementType(boundary::Boundary{T}, side::SIDE, index::Int)::TYPE where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
         throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
@@ -74,6 +164,20 @@ function getBoundaryElementType(boundary::Boundary{T}, side::SIDE, index::Int)::
     getType(boundary.boundaries[Int(side)][index])
 end
 
+"""
+    getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int)::T where {T}
+
+Retrieves the value of the boundary element at the specified side and index.
+Applicable if the boundary type is CONSTANT.
+
+# Arguments
+- `boundary::Boundary{T}`: The boundary.
+- `side::SIDE`: The side of the boundary.
+- `index::Int`: The index of the boundary element.
+
+# Returns
+The value of the boundary element.
+"""
 function getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int)::T where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
         throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
@@ -85,6 +189,18 @@ function getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int):
     getValue(boundary.boundaries[Int(side)][index])
 end
 
+"""
+    getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
+
+Retrieves the boundary elements at the specified side.
+
+# Arguments
+- `boundary::Boundary{T}`: The boundary.
+- `side::SIDE`: The side of the boundary.
+
+# Returns
+The boundary elements at the specified side.
+"""
 function getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
         throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
@@ -93,7 +209,19 @@ function getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElem
     boundary.boundaries[Int(side)]
 end
 
-function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE) where {T}
+"""
+    setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
+
+Sets the boundary elements at the specified side to CLOSED.
+
+# Arguments
+- `boundary::Boundary{T}`: The boundary.
+- `side::SIDE`: The side of the boundary.
+
+# Returns
+The boundary elements at the specified side.
+"""
+function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
         throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
     end
@@ -104,7 +232,20 @@ function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE) where {T}
     boundary.boundaries[Int(side)] = [BoundaryElement{T}() for _ in 1:n]
 end
 
-function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::T) where {T}
+"""
+    setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, values::Vector{T})::Vector{BoundaryElement{T}} where {T}
+
+Sets the boundary elements at the specified side to CONSTANT with the specified values.
+
+# Arguments
+- `boundary::Boundary{T}`: The boundary.
+- `side::SIDE`: The side of the boundary.
+- `value::T`: The value of the boundary elements.
+
+# Returns
+The boundary elements at the specified side.
+"""
+function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::T)::Vector{BoundaryElement{T}} where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
         throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
     end
diff --git a/julia/TUG/src/Core/BTCS.jl b/julia/TUG/src/Core/BTCS.jl
index 4662a1b..26a058a 100644
--- a/julia/TUG/src/Core/BTCS.jl
+++ b/julia/TUG/src/Core/BTCS.jl
@@ -8,78 +8,51 @@ using Base.Threads
 using LinearAlgebra
 using SparseArrays
 
-function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T}
+function calcBoundaryCoeff(alpha_center::T, alpha_side::T, sx::T, bcType::TYPE)::Tuple{T,T} where {T}
     alpha = calcAlphaIntercell(alpha_center, alpha_side)
-    centerCoeff = 1 + sx * alpha
     sideCoeff = -sx * alpha
-    return (centerCoeff, sideCoeff)
-end
 
-function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where {T}
-    alpha = calcAlphaIntercell(alpha_center, alpha_side)
-    centerCoeff = 1 + sx * (alpha + 2 * alpha_center)
-    sideCoeff = -sx * alpha
-    return (centerCoeff, sideCoeff)
-end
-
-# creates coefficient matrix for next time step from alphas in x-direction
-function createCoeffMatrix(alpha::Matrix{T}, alpha_left::Vector{T}, alpha_right::Vector{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, numCols::Int, rowIndex::Int, sx::T)::Tridiagonal{T} where {T}
-    # Precompute boundary condition type check for efficiency
-    bcTypeLeft = getType(bcLeft[rowIndex])
-
-    # Determine left side boundary coefficients based on boundary condition
-    centerCoeffTop, rightCoeffTop = if bcTypeLeft == CONSTANT
-        calcBoundaryCoeffConstant(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
-    elseif bcTypeLeft == CLOSED
-        calcBoundaryCoeffClosed(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
+    if bcType == CONSTANT
+        centerCoeff = 1 + sx * (alpha + 2 * alpha_center)
+    elseif bcType == CLOSED
+        centerCoeff = 1 + sx * alpha
     else
-        error("Undefined Boundary Condition Type on Left!")
+        error("Undefined Boundary Condition Type!")
     end
 
-    # Precompute boundary condition type check for efficiency
-    bcTypeRight = getType(bcRight[rowIndex])
+    return (centerCoeff, sideCoeff)
+end
+
+function createCoeffMatrix(alpha::Matrix{T}, alpha_left::Vector{T}, alpha_right::Vector{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, numCols::Int, rowIndex::Int, sx::T)::Tridiagonal{T} where {T}
+    # Determine left side boundary coefficients based on boundary condition
+    centerCoeffTop, rightCoeffTop = calcBoundaryCoeff(alpha[rowIndex, 1], alpha[rowIndex, 2], sx, getType(bcLeft[rowIndex]))
 
     # Determine right side boundary coefficients based on boundary condition
-    centerCoeffBottom, leftCoeffBottom = if bcTypeRight == CONSTANT
-        calcBoundaryCoeffConstant(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
-    elseif bcTypeRight == CLOSED
-        calcBoundaryCoeffClosed(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
-    else
-        error("Undefined Boundary Condition Type on Right!")
-    end
+    centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeff(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx, getType(bcRight[rowIndex]))
 
     dl = [-sx .* alpha_left; leftCoeffBottom]
     d = [centerCoeffTop; 1 .+ sx .* (alpha_right + alpha_left); centerCoeffBottom]
     du = [rightCoeffTop; -sx .* alpha_right]
-    alpha_diagonal = Tridiagonal(dl, d, du)
-
-    return alpha_diagonal
+    return Tridiagonal(dl, d, du)
 end
 
 
-function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
+function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T)::T where {T}
     alpha = calcAlphaIntercell(alpha_center, alpha_neighbor)
-
-    (sy * alpha + (1 - sy * alpha)) * conc_center
+    return (sy * alpha + (1 - sy * alpha)) * conc_center
 end
 
-function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
+function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::T)::T where {T}
     alpha_center_neighbor = calcAlphaIntercell(alpha_center, alpha_neighbor)
     alpha_center_center = alpha_center == alpha_neighbor ? alpha_center_neighbor : calcAlphaIntercell(alpha_center, alpha_center)
 
-    sy * alpha_center_neighbor * conc_center +
-    (1 - sy * (alpha_center_center + 2 * alpha_center)) * conc_center +
-    sy * alpha_center * conc_bc
+    return sy * alpha_center_neighbor * conc_center +
+           (1 - sy * (alpha_center_center + 2 * alpha_center)) * conc_center +
+           sy * alpha_center * conc_bc
 end
 
-function calcExplicitConcentrationsBoundaryClosed(conc_center::Vector{T}, alpha::Vector{T}, sy::T) where {T}
-    (sy .* alpha .+ (1 .- sy .* alpha)) .* conc_center
-end
-
-function calcExplicitConcentrationsBoundaryConstant(conc_center::Vector{T}, conc_bc::Vector{T}, alpha_center::Vector{T}, alpha_neighbor::Vector{T}, sy::T) where {T}
-    sy .* alpha_neighbor .* conc_center[rowIndex, :] +
-    (1 .- sy .* (alpha_center .+ 2 .* alphaY[rowIndex, :])) .* conc_center[rowIndex, :] +
-    sy .* alphaY[rowIndex, :] .* conc_bc
+function calcExplicitConcentrationsBoundaryClosed(conc_center::Vector{T}, alpha::Vector{T}, sy::T)::T where {T}
+    return (sy .* alpha .+ (1 .- sy .* alpha)) .* conc_center
 end
 
 # creates a solution vector for next time step from the current state of concentrations
@@ -87,17 +60,6 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
     numRows = size(concentrations, 1)
     length = size(sv, 1)
 
-    # Inner rows
-    if rowIndex > 1 && rowIndex < numRows
-        @inbounds for i = 1:length
-            alpha_here_below = calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i])
-            alpha_here_above = alphaY[rowIndex+1, i] == alphaY[rowIndex-1, i] ? alpha_here_below : calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) # calcAlphaIntercell is symmetric, so we can use it for both directions
-            sv[i] = sy * alpha_here_below * concentrations[rowIndex+1, i] +
-                    (1 - sy * (alpha_here_below + alpha_here_above)) * concentrations[rowIndex, i] +
-                    sy * alpha_here_above * concentrations[rowIndex-1, i]
-        end
-    end
-
     # First row
     if rowIndex == 1
         @inbounds for i = 1:length
@@ -111,6 +73,17 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
         end
     end
 
+    # Inner rows
+    if rowIndex > 1 && rowIndex < numRows
+        @inbounds for i = 1:length
+            alpha_here_below = calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i])
+            alpha_here_above = alphaY[rowIndex+1, i] == alphaY[rowIndex-1, i] ? alpha_here_below : calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) # calcAlphaIntercell is symmetric, so we can use it for both directions
+            sv[i] = sy * alpha_here_below * concentrations[rowIndex+1, i] +
+                    (1 - sy * (alpha_here_below + alpha_here_above)) * concentrations[rowIndex, i] +
+                    sy * alpha_here_above * concentrations[rowIndex-1, i]
+        end
+    end
+
     # Last row
     if rowIndex == numRows
         @inbounds for i = 1:length
@@ -135,9 +108,6 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
     end
 end
 
-
-
-# BTCS solution for 1D grid
 function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_right::Matrix{T}, timestep::T) where {T}
     sx = timestep / (getDeltaCol(grid) * getDeltaCol(grid))
 
@@ -161,7 +131,6 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_ri
     concentrations[1, :] = A \ b
 end
 
-# BTCS solution for 2D grid
 function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_right::Matrix{T}, alphaY_t_left::Matrix{T}, alphaY_t_right::Matrix{T}, timestep::T) where {T}
     rows = getRows(grid)
     cols = getCols(grid)
@@ -183,7 +152,6 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
     bcTop = getBoundarySide(bc, TOP)
     bcBottom = getBoundarySide(bc, BOTTOM)
 
-
     localBs = [zeros(T, cols) for _ in 1:Threads.nthreads()]
     Threads.@threads for i = 1:rows
         localB = localBs[Threads.threadid()]
@@ -193,11 +161,9 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
         concentrations_intermediate[i, :] = A \ localB
     end
 
-
     concentrations_intermediate = copy(transpose(concentrations_intermediate))
     concentrations_t = fetch(concentrations_t_task)
 
-
     # Swap alphas, boundary conditions and sx/sy for column-wise calculation
     localBs = [zeros(T, rows) for _ in 1:Threads.nthreads()]
     Threads.@threads for i = 1:cols
@@ -253,5 +219,4 @@ function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, s
     else
         error("BTCS only implemented for 1D and 2D grids")
     end
-
 end
diff --git a/julia/TUG/src/Core/FTCS.jl b/julia/TUG/src/Core/FTCS.jl
index 5daab41..26aede4 100644
--- a/julia/TUG/src/Core/FTCS.jl
+++ b/julia/TUG/src/Core/FTCS.jl
@@ -89,8 +89,6 @@ function calcVerticalChangeBottomBoundary(boundaryType::TYPE, alphaY_current::T,
     end
 end
 
-
-
 function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
     colMax = getCols(grid)
     sx = timestep / (getDeltaCol(grid)^2)
@@ -122,7 +120,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
     setConcentrations!(grid, concentrations_t1)
 end
 
-
 function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
     rowMax = getRows(grid)
     colMax = getCols(grid)
diff --git a/julia/TUG/src/Core/Utils.jl b/julia/TUG/src/Core/Utils.jl
index 36213a4..5ac6cc5 100644
--- a/julia/TUG/src/Core/Utils.jl
+++ b/julia/TUG/src/Core/Utils.jl
@@ -2,10 +2,6 @@ function calcAlphaIntercell(alpha1::T, alpha2::T) where {T}
     2 / ((1 / alpha1) + (1 / alpha2))
 end
 
-function calcAlphaIntercell(alpha1::Vector{T}, alpha2::Vector{T}) where {T}
-    2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
-end
-
 function calcAlphaIntercell(alpha1::Matrix{T}, alpha2::Matrix{T}) where {T}
     2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
 end
diff --git a/julia/TUG/src/DynamicSimulation.jl b/julia/TUG/src/DynamicSimulation.jl
index 9a55661..756566d 100644
--- a/julia/TUG/src/DynamicSimulation.jl
+++ b/julia/TUG/src/DynamicSimulation.jl
@@ -1,5 +1,22 @@
 using Distributed
 
+"""
+    struct DynamicSimulation{T} <: AbstractSimulation{T}
+
+Represents a dynamic simulation environment with multiple grid states. 
+Allows the manipulation and running of simulations on different grid states in parallel.
+
+# Fields
+- `grid::Grid{T}`: The initial grid state.
+- `grids::Vector{Grid{T}}`: A vector of grid states for the simulation.
+- `bc::Boundary{T}`: The boundary conditions for the simulation.
+- `approach::APPROACH`: The numerical approach (BTCS or FTCS).
+- `iterations::Int`: Number of iterations to run each simulation.
+- `timestep::T`: The timestep for each iteration.
+
+# Constructor
+- `DynamicSimulation(grid, bc, approach, timestep)` creates a new dynamic simulation with specified parameters.
+"""
 struct DynamicSimulation{T} <: AbstractSimulation{T}
     grid::Grid{T}
     grids::Vector{Grid{T}}
@@ -9,34 +26,79 @@ struct DynamicSimulation{T} <: AbstractSimulation{T}
     iterations::Int
     timestep::T
 
-    # Constructor
-    function DynamicSimulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH, timestep::T) where {T}
+    function DynamicSimulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH, timestep::T)::DynamicSimulation{T} where {T}
         timestep, iterations = adjustTimestep(grid, approach, timestep, 1, false)
         new{T}(grid, Vector{Grid{T}}(), bc, approach, iterations, timestep)
     end
 end
 
+"""
+    createGrid(simulation::DynamicSimulation{T})::Int where {T}
+
+Creates a new grid state based on the initial state and adds it to the simulation's grid vector.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+
+# Returns
+The index of the newly created grid in the simulation's grid vector.
+"""
 function createGrid(simulation::DynamicSimulation{T})::Int where {T}
     new_grid = clone(simulation.grid)
     push!(simulation.grids, new_grid)
     return length(simulation.grids)
 end
 
-function next(simulation::DynamicSimulation{T}) where {T}
+"""
+    next(simulation::DynamicSimulation{T})::Nothing where {T}
+
+Runs the simulation for all grid states in parallel.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+
+# Returns
+Nothing, but updates all grids in the simulation according to the defined approach and iterations.
+"""
+function next(simulation::DynamicSimulation{T})::Nothing where {T}
     pmap(grid -> runSimulationForGrid(simulation, grid), simulation.grids)
 end
 
-function printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int) where {T}
+"""
+    printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
+
+Prints the concentrations of a specific grid state to the console.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
+
+# Returns
+Nothing, but outputs the grid's concentration matrix to the console.
+"""
+function printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
     writeConcentrationsToCLI(simulation, simulation.grids[gridIndex])
 end
 
-function printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int) where {T}
+"""
+    printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
+
+Writes the concentrations of a specific grid state to a CSV file.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
+
+# Returns
+Nothing, but writes the grid's concentration matrix to a CSV file.
+"""
+function printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
     file = createCSVfile(simulation, simulation.grids[gridIndex])
     writeConcentrationsToCSV(file, simulation, simulation.grids[gridIndex])
     close(file)
 end
 
-function runSimulationForGrid(simulation::DynamicSimulation{T}, grid::Grid{T}) where {T}
+function runSimulationForGrid(simulation::DynamicSimulation{T}, grid::Grid{T})::Nothing where {T}
     if simulation.approach == BTCS
         runBTCS(grid, simulation.bc, simulation.timestep, simulation.iterations, () -> nothing)
     elseif simulation.approach == FTCS
@@ -46,18 +108,69 @@ function runSimulationForGrid(simulation::DynamicSimulation{T}, grid::Grid{T}) w
     end
 end
 
+"""
+    getConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Matrix{T} where {T}
+
+Retrieves the concentrations of a specific grid state.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
+
+# Returns
+The concentration matrix of the specified grid.
+"""
 function getConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Matrix{T} where {T}
     getConcentrations(simulation.grids[gridIndex])
 end
 
-function setConcentrations!(simulation::DynamicSimulation{T}, gridIndex::Int, concentrations::Matrix{T}) where {T}
+"""
+    setConcentrations!(simulation::DynamicSimulation{T}, gridIndex::Int, concentrations::Matrix{T})::Nothing where {T}
+
+Sets the concentrations of a specific grid state.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
+- `concentrations::Matrix{T}`: The new concentration matrix.
+
+# Returns
+Nothing, but updates the concentration matrix of the specified grid.
+"""
+function setConcentrations!(simulation::DynamicSimulation{T}, gridIndex::Int, concentrations::Matrix{T})::Nothing where {T}
     setConcentrations!(simulation.grids[gridIndex], concentrations)
 end
 
-function setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T}) where {T}
+"""
+    setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T})::Nothing where {T}
+
+Sets the alphaX matrix of a specific grid state.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
+- `alphaX::Matrix{T}`: The new alphaX matrix.
+
+# Returns
+Nothing, but updates the alphaX matrix of the specified grid.
+"""
+function setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T})::Nothing where {T}
     setAlphaX!(simulation.grids[gridIndex], alphaX)
 end
 
-function setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T}) where {T}
+"""
+    setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T})::Nothing where {T}
+
+Sets the alphaY matrix of a specific grid state.
+
+# Arguments
+- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
+- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
+- `alphaY::Matrix{T}`: The new alphaY matrix.
+
+# Returns
+Nothing, but updates the alphaY matrix of the specified grid.
+"""
+function setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T})::Nothing where {T}
     setAlphaY!(simulation.grids[gridIndex], alphaY)
 end
diff --git a/julia/TUG/src/Grid.jl b/julia/TUG/src/Grid.jl
index ae47451..321051f 100644
--- a/julia/TUG/src/Grid.jl
+++ b/julia/TUG/src/Grid.jl
@@ -5,6 +5,28 @@
 
 using LinearAlgebra
 
+"""
+    struct Grid{T}
+
+Represents a computational grid with concentration values and alpha coefficients.
+It can be either 1D or 2D and includes dimensions, domain sizes, delta values, and matrices for concentrations and alpha coefficients.
+
+# Fields
+- `cols::Int`: Number of columns in the grid.
+- `rows::Int`: Number of rows in the grid.
+- `dim::Int`: Dimension of the grid (1 for 1D, 2 for 2D).
+- `domainCol::T`, `domainRow::T`: Size of the grid's domain in column and row direction.
+- `deltaCol::T`, `deltaRow::T`: Delta values for columns and rows.
+- `concentrations::Ref{Matrix{T}}`: Reference to the matrix holding concentration values.
+- `alphaX::Ref{Matrix{T}}`: Reference to the matrix of alpha coefficients in the X direction.
+- `alphaY::Union{Ref{Matrix{T}},Nothing}`: Reference to the matrix of alpha coefficients in the Y direction (if applicable).
+- `alphaX_t`, `alphaY_t`: Transposed matrices of alpha coefficients.
+
+# Constructors
+- `Grid(length, alphaX)` creates a 1D grid with the given length and alphaX matrix.
+- `Grid(rows, cols, alphaX, alphaY)` creates a 2D grid with the given rows, columns, and alphaX and alphaY matrices.
+- `Grid(rows, cols, dim, domainCol, domainRow, deltaCol, deltaRow, concentrations, alphaX, alphaY, alphaX_t, alphaY_t)` creates a grid with the given parameters.
+"""
 struct Grid{T}
     cols::Int
     rows::Int
@@ -19,8 +41,7 @@ struct Grid{T}
     alphaX_t::Union{Ref{Matrix{T}},Nothing}
     alphaY_t::Union{Ref{Matrix{T}},Nothing}
 
-    # Constructor for 1D-Grid
-    function Grid{T}(length::Int, alphaX::Matrix{T}) where {T}
+    function Grid{T}(length::Int, alphaX::Matrix{T})::Grid{T} where {T}
         if length <= 3
             throw(ArgumentError("Given grid length too small. Must be greater than 3."))
         end
@@ -33,8 +54,7 @@ struct Grid{T}
         new{T}(length, 1, 1, T(length), 0, T(1), 0, Ref(fill(T(0), 1, length)), alphaX, nothing, alphaX_t, nothing)
     end
 
-    # Constructor for 2D-Grid
-    function Grid{T}(rows::Int, cols::Int, alphaX::Matrix{T}, alphaY::Matrix{T}) where {T}
+    function Grid{T}(rows::Int, cols::Int, alphaX::Matrix{T}, alphaY::Matrix{T})::Grid{T} where {T}
         if rows <= 3 || cols <= 3
             throw(ArgumentError("Given grid dimensions too small. Must each be greater than 3."))
         end
@@ -42,18 +62,29 @@ struct Grid{T}
             throw(ArgumentError("Given matrices of alpha coefficients mismatch with Grid dimensions!"))
         end
 
-        # Precompute alphaX_t and alphaY_t
         alphaX_t = alphaX'
         alphaY_t = alphaY'
 
         new{T}(cols, rows, 2, T(cols), T(rows), T(1), T(1), Ref(fill(T(0), rows, cols)), alphaX, alphaY, alphaX_t, alphaY_t)
     end
 
-    function Grid{T}(rows::Int, cols::Int, dim::Int, domainCol::T, domainRow::T, deltaCol::T, deltaRow::T, concentrations::Ref{Matrix{T}}, alphaX::Ref{Matrix{T}}, alphaY::Union{Ref{Matrix{T}},Nothing}, alphaX_t::Union{Ref{Matrix{T}},Nothing}, alphaY_t::Union{Ref{Matrix{T}},Nothing}) where {T}
+    function Grid{T}(rows::Int, cols::Int, dim::Int, domainCol::T, domainRow::T, deltaCol::T, deltaRow::T, concentrations::Ref{Matrix{T}}, alphaX::Ref{Matrix{T}}, alphaY::Union{Ref{Matrix{T}},Nothing}, alphaX_t::Union{Ref{Matrix{T}},Nothing}, alphaY_t::Union{Ref{Matrix{T}},Nothing})::Grid{T} where {T}
         new{T}(cols, rows, dim, domainCol, domainRow, deltaCol, deltaRow, concentrations, alphaX, alphaY, alphaX_t, alphaY_t)
     end
 end
 
+"""
+    clone(grid::Grid{T})::Grid{T} where {T}
+
+Creates a deep copy of the specified grid. The new grid will have identical properties 
+but will be a distinct object in memory.
+
+# Arguments
+- `grid::Grid{T}`: The grid to be cloned.
+
+# Returns
+A new `Grid{T}` object that is a deep copy of the input grid.
+"""
 function clone(grid::Grid{T})::Grid{T} where {T}
     if grid.dim == 1
         return Grid{T}(1, grid.cols, grid.dim, grid.domainCol, grid.domainRow, grid.deltaCol, grid.deltaRow, Ref(copy(grid.concentrations[])), Ref(copy(grid.alphaX[])), nothing, Ref(copy(grid.alphaX_t[])), nothing)
@@ -61,10 +92,33 @@ function clone(grid::Grid{T})::Grid{T} where {T}
     Grid{T}(grid.rows, grid.cols, grid.dim, grid.domainCol, grid.domainRow, grid.deltaCol, grid.deltaRow, Ref(copy(grid.concentrations[])), Ref(copy(grid.alphaX[])), Ref(copy(grid.alphaY[])), Ref(copy(grid.alphaX_t[])), Ref(copy(grid.alphaY_t[])))
 end
 
+"""
+    getAlphaX(grid::Grid{T})::Matrix{T} where {T}
+
+Retrieves the alpha coefficients matrix in the X direction from the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the alphaX matrix.
+
+# Returns
+The `alphaX` matrix of the grid.
+"""
 function getAlphaX(grid::Grid{T})::Matrix{T} where {T}
     grid.alphaX[]
 end
 
+"""
+    getAlphaY(grid::Grid{T})::Matrix{T} where {T}
+
+Retrieves the alpha coefficients matrix in the Y direction from the specified grid.
+Raises an error if the grid is 1D (no alphaY matrix).
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the alphaY matrix.
+
+# Returns
+The `alphaY` matrix of the grid, if applicable.
+"""
 function getAlphaY(grid::Grid{T})::Matrix{T} where {T}
     if grid.dim == 1
         error("Grid is 1D, so there is no alphaY matrix!")
@@ -73,10 +127,33 @@ function getAlphaY(grid::Grid{T})::Matrix{T} where {T}
     grid.alphaY[]
 end
 
+"""
+    getAlphaX_t(grid::Grid{T})::Matrix{T} where {T}
+
+Retrieves the transposed alpha coefficients matrix in the X direction from the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the alphaX_t matrix.
+
+# Returns
+The transposed `alphaX_t` matrix of the grid.
+"""
 function getAlphaX_t(grid::Grid{T})::Matrix{T} where {T}
     grid.alphaX_t[]
 end
 
+"""
+    getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
+
+Retrieves the transposed alpha coefficients matrix in the Y direction from the specified grid.
+Raises an error if the grid is 1D (no alphaY_t matrix).
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the alphaY_t matrix.
+
+# Returns
+The transposed `alphaY_t` matrix of the grid, if applicable.
+"""
 function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
     if grid.dim == 1
         error("Grid is 1D, so there is no alphaY_t matrix!")
@@ -85,40 +162,133 @@ function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
     grid.alphaY_t[]
 end
 
+"""
+    getCols(grid::Grid{T})::Int where {T}
+
+Retrieves the number of columns in the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the number of columns.
+
+# Returns
+The number of columns in the grid.
+"""
 function getCols(grid::Grid{T})::Int where {T}
     grid.cols
 end
 
+"""
+    getConcentrations(grid::Grid{T})::Matrix{T} where {T}
+
+Retrieves the concentration matrix from the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the concentration matrix.
+
+# Returns
+The concentration matrix of the grid.
+"""
 function getConcentrations(grid::Grid{T})::Matrix{T} where {T}
     grid.concentrations[]
 end
 
+"""
+    getDeltaCol(grid::Grid{T})::T where {T}
+
+Retrieves the delta value for columns in the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the delta column value.
+
+# Returns
+The delta value for columns in the grid.
+"""
 function getDeltaCol(grid::Grid{T})::T where {T}
     grid.deltaCol
 end
 
+"""
+    getDeltaRow(grid::Grid{T})::T where {T}
+
+Retrieves the delta value for rows in the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the delta row value.
+
+# Returns
+The delta value for rows in the grid.
+"""
 function getDeltaRow(grid::Grid{T})::T where {T}
     grid.deltaRow
 end
 
+"""
+    getDim(grid::Grid{T})::Int where {T}
+
+Retrieves the dimension (1D or 2D) of the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the dimension.
+
+# Returns
+The dimension of the grid (1 for 1D, 2 for 2D).
+"""
 function getDim(grid::Grid{T})::Int where {T}
     grid.dim
 end
 
+"""
+    getRows(grid::Grid{T})::Int where {T}
+
+Retrieves the number of rows in the specified grid.
+
+# Arguments
+- `grid::Grid{T}`: The grid from which to retrieve the number of rows.
+
+# Returns
+The number of rows in the grid.
+"""
 function getRows(grid::Grid{T})::Int where {T}
     grid.rows
 end
 
-function setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T}) where {T}
+"""
+    setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T})::Nothing where {T}
+
+Sets the alpha coefficients matrix in the X direction for the specified grid.
+Throws an error if the dimensions of the new matrix don't match the grid's dimensions.
+
+# Arguments
+- `grid::Grid{T}`: The grid for which to set the new alphaX matrix.
+- `new_alphaX::Matrix{T}`: The new alpha coefficients matrix in the X direction.
+
+# Returns
+Nothing, but modifies the alphaX matrix of the grid.
+"""
+function setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T})::Nothing where {T}
     if size(new_alphaX) != size(grid.alphaX[])
         throw(ArgumentError("Given matrix of alpha coefficients mismatch with Grid dimensions!"))
     end
 
     grid.alphaX[] = new_alphaX
     grid.alphaX_t[] = new_alphaX'
+    return
 end
 
-function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T}) where {T}
+"""
+    setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T})::Nothing where {T}
+
+Sets the alpha coefficients matrix in the Y direction for the specified grid.
+Throws an error if the grid is 1D or if the dimensions of the new matrix don't match the grid's dimensions.
+
+# Arguments
+- `grid::Grid{T}`: The grid for which to set the new alphaY matrix.
+- `new_alphaY::Matrix{T}`: The new alpha coefficients matrix in the Y direction.
+
+# Returns
+Nothing, but modifies the alphaY matrix of the grid.
+"""
+function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T})::Nothing where {T}
     if grid.dim == 1
         error("Grid is 1D, so there is no alphaY matrix!")
     end
@@ -128,12 +298,27 @@ function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T}) where {T}
 
     grid.alphaY[] = new_alphaY
     grid.alphaY_t[] = new_alphaY'
+    return
 end
 
-function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where {T}
+"""
+    setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T})::Nothing where {T}
+
+Sets the concentration matrix for the specified grid.
+Throws an error if the dimensions of the new matrix don't match the grid's dimensions.
+
+# Arguments
+- `grid::Grid{T}`: The grid for which to set the new concentrations matrix.
+- `new_concentrations::Matrix{T}`: The new concentrations matrix.
+
+# Returns
+Nothing, but modifies the concentration matrix of the grid.
+"""
+function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T})::Nothing where {T}
     if size(new_concentrations) != size(grid.concentrations[])
         throw(ArgumentError("Given matrix of concentrations mismatch with Grid dimensions!"))
     end
 
     grid.concentrations[] = new_concentrations
+    return
 end
diff --git a/julia/TUG/src/Simulation.jl b/julia/TUG/src/Simulation.jl
index 82e3f34..6e27542 100644
--- a/julia/TUG/src/Simulation.jl
+++ b/julia/TUG/src/Simulation.jl
@@ -4,6 +4,25 @@
 # options. Simulation object also holds a predefined Grid and Boundary object.
 # Translated from C++'s Simulation.hpp.
 
+"""
+    struct Simulation{T} <: AbstractSimulation{T}
+
+Represents a simulation environment with a specified grid and boundary conditions. 
+Holds all necessary information for a simulation run including timestep, number 
+of iterations, output options, and simulation approach.
+
+# Fields
+- `grid::Grid{T}`: Grid object defining the simulation space.
+- `bc::Boundary{T}`: Boundary conditions for the simulation.
+- `approach::APPROACH`: Numerical approach for the simulation (e.g., BTCS, FTCS).
+- `iterations::Int`: Number of iterations the simulation will run.
+- `timestep::T`: Time step for each simulation iteration.
+- `consoleOutput::CONSOLE_OUTPUT`: Option for console output level.
+- `csvOutput::CSV_OUTPUT`: Option for CSV file output level.
+
+# Constructor
+- `Simulation(grid, bc, approach, iterations, timestep, consoleOutput, csvOutput)` creates a new simulation with specified parameters.
+"""
 struct Simulation{T} <: AbstractSimulation{T}
     grid::Grid{T}
     bc::Boundary{T}
@@ -15,14 +34,25 @@ struct Simulation{T} <: AbstractSimulation{T}
     consoleOutput::CONSOLE_OUTPUT
     csvOutput::CSV_OUTPUT
 
-    # Constructor
     function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS, iterations::Int=1, timestep::T=0.1,
-        consoleOutput::CONSOLE_OUTPUT=CONSOLE_OUTPUT_OFF, csvOutput::CSV_OUTPUT=CSV_OUTPUT_OFF) where {T}
+        consoleOutput::CONSOLE_OUTPUT=CONSOLE_OUTPUT_OFF, csvOutput::CSV_OUTPUT=CSV_OUTPUT_OFF)::Simulation{T} where {T}
         new{T}(grid, bc, approach, iterations, timestep, consoleOutput, csvOutput)
     end
 end
 
-function run(simulation::Simulation{T}) where {T}
+"""
+    run(simulation::Simulation{T})::Nothing where {T}
+
+Executes the simulation based on the defined parameters within the `simulation` object.
+Depending on the output settings, it may write data to console and/or a CSV file.
+
+# Arguments
+- `simulation::Simulation{T}`: The simulation object to be run.
+
+# Returns
+Nothing, but executes the simulation and handles outputs.
+"""
+function run(simulation::Simulation{T})::Nothing where {T}
     file = nothing
     try
         if simulation.csvOutput > CSV_OUTPUT_OFF
@@ -54,7 +84,6 @@ function run(simulation::Simulation{T}) where {T}
         if simulation.csvOutput >= CSV_OUTPUT_ON
             writeConcentrationsToCSV(file, simulation)
         end
-
     finally
         if file !== nothing
             close(file)
@@ -62,18 +91,66 @@ function run(simulation::Simulation{T}) where {T}
     end
 end
 
+"""
+    setIterations(simulation::Simulation{T}, iterations::Int)::Simulation{T} where {T}
+
+Sets the number of iterations for the given simulation.
+
+# Arguments
+- `simulation::Simulation{T}`: The simulation object.
+- `iterations::Int`: The new number of iterations to be set.
+
+# Returns
+A new `Simulation` object with updated iterations.
+"""
 function setIterations(simulation::Simulation{T}, iterations::Int)::Simulation{T} where {T}
     return Simulation(simulation.grid, simulation.bc, simulation.approach, iterations, simulation.timestep, simulation.consoleOutput, simulation.csvOutput)
 end
 
+"""
+    setOutputConsole(simulation::Simulation{T}, consoleOutput::CONSOLE_OUTPUT)::Simulation{T} where {T}
+
+Sets the console output level for the simulation.
+
+# Arguments
+- `simulation::Simulation{T}`: The simulation object.
+- `consoleOutput::CONSOLE_OUTPUT`: The new console output level.
+
+# Returns
+A new `Simulation` object with updated console output setting.
+"""
 function setOutputConsole(simulation::Simulation{T}, consoleOutput::CONSOLE_OUTPUT)::Simulation{T} where {T}
     return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, consoleOutput, simulation.csvOutput)
 end
 
+"""
+    setOutputCSV(simulation::Simulation{T}, csvOutput::CSV_OUTPUT)::Simulation{T} where {T}
+
+Sets the CSV output level for the simulation.
+
+# Arguments
+- `simulation::Simulation{T}`: The simulation object.
+- `csvOutput::CSV_OUTPUT`: The new CSV output level.
+
+# Returns
+A new `Simulation` object with updated CSV output setting.
+"""
 function setOutputCSV(simulation::Simulation{T}, csvOutput::CSV_OUTPUT)::Simulation{T} where {T}
     return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, simulation.consoleOutput, csvOutput)
 end
 
+"""
+    setTimestep(simulation::Simulation{T}, timestep::T)::Simulation{T} where {T}
+
+Sets the timestep for the simulation.
+
+# Arguments
+- `simulation::Simulation{T}`: The simulation object.
+- `timestep::T`: The new timestep to be set.
+
+# Returns
+A new `Simulation` object with updated timestep.
+"""
 function setTimestep(simulation::Simulation{T}, timestep::T)::Simulation{T} where {T}
     return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, timestep, simulation.consoleOutput, simulation.csvOutput)
 end
diff --git a/julia/TUG/src/TUG.jl b/julia/TUG/src/TUG.jl
index 723e230..fa2afd4 100644
--- a/julia/TUG/src/TUG.jl
+++ b/julia/TUG/src/TUG.jl
@@ -3,7 +3,7 @@ module TUG
 include("Grid.jl")
 
 export Grid
-export getAlphaX, getAlphaY, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!, getDomainCol, getDomainRow, getDeltaCol, getDeltaRow, getDim
+export clone, getAlphaX, getAlphaY, getAlphaX_t, getAlphaY_t, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!, getDomainCol, getDomainRow, getDeltaCol, getDeltaRow, getDim
 
 include("Boundary.jl")
 
@@ -12,7 +12,7 @@ export getType, getValue, setBoundarySideClosed!, setBoundarySideConstant!, getB
 
 include("AbstractSimulation.jl")
 
-export AbstractSimulation, APPROACH, CONSOLE_OUTPUT, CSV_OUTPUT, BTCS, FTCS, CONSOLE_OUTPUT_OFF, CONSOLE_OUTPUT_ON, CONSOLE_OUTPUT_VERBOSE, CSV_OUTPUT_OFF, CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE, CSV_OUTPUT_XTREME
+export APPROACH, CONSOLE_OUTPUT, CSV_OUTPUT, BTCS, FTCS, CONSOLE_OUTPUT_OFF, CONSOLE_OUTPUT_ON, CONSOLE_OUTPUT_VERBOSE, CSV_OUTPUT_OFF, CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE, CSV_OUTPUT_XTREME
 
 include("Simulation.jl")
 
@@ -22,7 +22,7 @@ export run, setTimestep, setIterations, setOutputConsole, setOutputCSV
 include("DynamicSimulation.jl")
 
 export DynamicSimulation
-export createGrid, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!, next, printConcentrationsCSV, printConcentrations
+export createGrid, next, printConcentrations, printConcentrationsCSV, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!
 
 include("Core/Utils.jl")
 include("Core/BTCS.jl")