diff --git a/julia/tug/Boundary.jl b/julia/tug/Boundary.jl
index 455f7e8..cb6cb60 100644
--- a/julia/tug/Boundary.jl
+++ b/julia/tug/Boundary.jl
@@ -1,5 +1,7 @@
 # Boundary.jl
-# Julia implementation of Boundary types, translating from C++'s Boundary.hpp
+# API of Boundary class, that holds all information for each boundary
+# condition at the edges of the diffusion grid. 
+# Translated from C++'s Boundary.hpp.
 
 @enum TYPE CLOSED CONSTANT
 @enum SIDE LEFT = 1 RIGHT = 2 TOP = 3 BOTTOM = 4
@@ -17,12 +19,6 @@ struct BoundaryElement{T}
     BoundaryElement{T}(value::T) where {T} = new{T}(CONSTANT, value)
 end
 
-
-# Helper functions for setting type and value
-function setType!(be::BoundaryElement{T}, type::Symbol) where {T}
-    be.type = type
-end
-
 function getType(be::BoundaryElement{T})::TYPE where {T}
     be.type
 end
@@ -31,6 +27,10 @@ function getValue(be::BoundaryElement{T})::T where {T}
     be.value
 end
 
+function setType!(be::BoundaryElement{T}, type::Symbol) where {T}
+    be.type = type
+end
+
 function setValue!(be::BoundaryElement{T}, value::T) where {T}
     if value < 0
         throw(ArgumentError("No negative concentration allowed."))
@@ -38,9 +38,11 @@ function setValue!(be::BoundaryElement{T}, value::T) where {T}
     if be.type == BC_TYPE_CLOSED
         throw(ArgumentError("No constant boundary concentrations can be set for closed boundaries. Please change type first."))
     end
+
     be.value = value
 end
 
+
 # Boundary class
 struct Boundary{T}
     dim::UInt8
@@ -71,6 +73,14 @@ struct Boundary{T}
     end
 end
 
+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."))
+    end
+
+    boundary.boundaries[Int(side)]
+end
+
 function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE) where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
         throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
@@ -92,11 +102,3 @@ function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::T) w
 
     boundary.boundaries[Int(side)] = [BoundaryElement{T}(value) for _ in 1:n]
 end
-
-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."))
-    end
-
-    boundary.boundaries[Int(side)]
-end
diff --git a/julia/tug/Core/BTCS.jl b/julia/tug/Core/BTCS.jl
index 648aa2a..cb51fd0 100644
--- a/julia/tug/Core/BTCS.jl
+++ b/julia/tug/Core/BTCS.jl
@@ -2,21 +2,29 @@
 # Implementation of heterogenous BTCS (backward time-centered space)
 # solution of diffusion equation in 1D and 2D space using the
 # alternating-direction implicit (ADI) method.
+# Translated from C++'s BTCS.hpp.
 
 using LinearAlgebra
 using SparseArrays
 using Base.Threads
-using CUDA
 
 include("../Boundary.jl")
 include("../Grid.jl")
 
 function calcAlphaIntercell(alpha1::T, alpha2::T) where {T}
-    return 2 / ((1 / alpha1) + (1 / alpha2))
+    2 / ((1 / alpha1) + (1 / alpha2))
 end
 
 function calcAlphaIntercell(alpha1::Matrix{T}, alpha2::Matrix{T}) where {T}
-    return 2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
+    2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
+end
+
+
+function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::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}
@@ -26,13 +34,6 @@ function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where
     return (centerCoeff, sideCoeff)
 end
 
-function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T}
-    alpha = calcAlphaIntercell(alpha_center, alpha_side)
-    centerCoeff = 1 + sx * alpha
-    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
@@ -70,18 +71,20 @@ end
 
 function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
     alpha = calcAlphaIntercell(alpha_center, alpha_neighbor)
-    sy * alpha * conc_center + (1 - sy * alpha) * conc_center
-end
 
+    (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}
     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
 end
 
+# creates a solution vector for next time step from the current state of concentrations inplace
 function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::Matrix{T}, alphaY::Matrix{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, bcTop::Vector{BoundaryElement{T}}, bcBottom::Vector{BoundaryElement{T}}, rowIndex::Int, sx::T, sy::T) where {T}
     numRows = size(concentrations, 1)
     length = size(sv, 1)
@@ -134,9 +137,10 @@ 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 / (grid.deltaCol * grid.deltaCol)
+    sx = timestep / (getDeltaCol(grid) * getDeltaCol(grid))
 
     alpha = getAlphaX(grid)
     bcLeft = getBoundarySide(bc, LEFT)
@@ -163,8 +167,8 @@ end
 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)
-    sx = timestep / (2 * grid.deltaCol * grid.deltaCol)
-    sy = timestep / (2 * grid.deltaRow * grid.deltaRow)
+    sx = timestep / (2 * getDeltaCol(grid) * getDeltaCol(grid))
+    sy = timestep / (2 * getDeltaRow(grid) * getDeltaRow(grid))
 
     alphaX = getAlphaX(grid)
     alphaY = getAlphaY(grid)
@@ -209,7 +213,7 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
 end
 
 function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, stepCallback::Function) where {T}
-    if grid.dim == 1
+    if getDim(grid) == 1
         length = getCols(grid)
         alpha = getAlphaX(grid)
 
@@ -224,7 +228,7 @@ function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, s
 
             stepCallback()
         end
-    elseif grid.dim == 2
+    elseif getDim(grid) == 2
         rows = getRows(grid)
         cols = getCols(grid)
         alphaX = getAlphaX(grid)
diff --git a/julia/tug/Grid.jl b/julia/tug/Grid.jl
index 85b4d50..884bf9c 100644
--- a/julia/tug/Grid.jl
+++ b/julia/tug/Grid.jl
@@ -1,5 +1,11 @@
+# Grid.jl
+# API of Grid class, that holds a matrix with concenctrations and 
+# a respective matrix/matrices of alpha coefficients.
+# Translated from C++'s Grid.hpp.
+
 using LinearAlgebra
 
+# Grid class
 struct Grid{T}
     cols::Int
     rows::Int
@@ -35,6 +41,7 @@ struct Grid{T}
             error("Given matrices of alpha coefficients mismatch with Grid dimensions!")
         end
 
+        # Precompute alphaX_t and alphaY_t
         alphaX_t = alphaX'
         alphaY_t = alphaY'
 
@@ -58,18 +65,30 @@ function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
     grid.alphaY_t
 end
 
+function getCols(grid::Grid{T})::Int where {T}
+    grid.cols
+end
+
 function getConcentrations(grid::Grid{T})::Matrix{T} where {T}
     grid.concentrations[]
 end
 
-function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where {T}
-    grid.concentrations[] = new_concentrations
+function getDeltaCol(grid::Grid{T})::T where {T}
+    grid.deltaCol
 end
 
-function getCols(grid::Grid{T})::Int where {T}
-    grid.cols
+function getDeltaRow(grid::Grid{T})::T where {T}
+    grid.deltaRow
+end
+
+function getDim(grid::Grid{T})::Int where {T}
+    grid.dim
 end
 
 function getRows(grid::Grid{T})::Int where {T}
     grid.rows
 end
+
+function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where {T}
+    grid.concentrations[] = new_concentrations
+end
diff --git a/julia/tug/Simulation.jl b/julia/tug/Simulation.jl
index 1278974..aa8a03b 100644
--- a/julia/tug/Simulation.jl
+++ b/julia/tug/Simulation.jl
@@ -1,3 +1,9 @@
+# Simulation.jl
+# API of Simulation class, that holds all information regarding a
+# specific simulation run like its timestep, number of iterations and output
+# options. Simulation object also holds a predefined Grid and Boundary object.
+# Translated from C++'s Simulation.hpp.
+
 using Printf
 
 include("Grid.jl")
@@ -5,16 +11,14 @@ include("Boundary.jl")
 include("Core/BTCS.jl")
 
 @enum APPROACH BTCS
-@enum SOLVER EIGEN_LU_SOLVER
 @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
 
-# Create the Simulation class
-struct Simulation{T,approach,solver}
+# Simulation class
+struct Simulation{T}
     grid::Grid{T}
     bc::Boundary{T}
     approach::APPROACH
-    solver::SOLVER
 
     iterations::Int
     timestep::T
@@ -22,36 +26,36 @@ struct Simulation{T,approach,solver}
     consoleOutput::CONSOLE_OUTPUT
     csvOutput::CSV_OUTPUT
 
-    function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS,
-        solver::SOLVER=EIGEN_LU_SOLVER, iterations::Int=1, timestep::T=0.1,
+    # 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}
-        new{T,APPROACH,SOLVER}(grid, bc, approach, solver, iterations, timestep, consoleOutput, csvOutput)
+        new{T}(grid, bc, approach, iterations, timestep, consoleOutput, csvOutput)
     end
-
 end
 
-function createCSVfile(simulation::Simulation{T,approach,solver})::IOStream where {T,approach,solver}
-    appendIdent = 0
-    approachString = (simulation.approach == BTCS) ? "BTCS" : "UNKNOWN"  # Add other approaches as needed
-    row = simulation.grid.rows
-    col = simulation.grid.cols
+function _createCSVfile(simulation::Simulation{T})::IOStream where {T}
+    approachString = string(simulation.approach)
+    rows = getRows(simulation.grid)
+    cols = getCols(simulation.grid)
     numIterations = simulation.iterations
-    filename = string(approachString, "_", row, "_", col, "_", numIterations, ".csv")
 
+    filename = string(approachString, "_", rows, "_", cols, "_", numIterations, ".csv")
+
+    appendIdent = 0
     while isfile(filename)
         appendIdent += 1
-        filename = string(approachString, "_", row, "_", col, "_", numIterations, "-", appendIdent, ".csv")
+        filename = string(approachString, "_", rows, "_", cols, "_", numIterations, "-", appendIdent, ".csv")
     end
 
     # Write boundary conditions if required
-    if simulation.csvOutput == CSV_OUTPUT_XTREME
+    if simulation.csvOutput >= CSV_OUTPUT_XTREME
         open(filename, "w") do file
-            writeBoundarySideValues(file, simulation.bc, LEFT)
-            writeBoundarySideValues(file, simulation.bc, RIGHT)
+            _writeBoundarySideValues(file, simulation.bc, LEFT)
+            _writeBoundarySideValues(file, simulation.bc, RIGHT)
 
-            if simulation.grid.dim == 2
-                writeBoundarySideValues(file, simulation.bc, TOP)
-                writeBoundarySideValues(file, simulation.bc, BOTTOM)
+            if getDim(simulation.grid) == 2
+                _writeBoundarySideValues(file, simulation.bc, TOP)
+                _writeBoundarySideValues(file, simulation.bc, BOTTOM)
             end
 
             write(file, "\n\n")
@@ -62,40 +66,40 @@ function createCSVfile(simulation::Simulation{T,approach,solver})::IOStream wher
     return file
 end
 
-function writeBoundarySideValues(file, bc::Boundary{T}, side) where {T}
+function _writeBoundarySideValues(file::IOStream, bc::Boundary{T}, side) where {T}
     values::Vector{BoundaryElement} = getBoundarySide(bc, side)
     formatted_values = join(map(getValue, values), " ")
     write(file, formatted_values, "\n")
 end
 
-function printConcentrationsCSV(simulation::Simulation{T,approach,solver}, file::IOStream) where {T,approach,solver}
-    concentrations = simulation.grid.concentrations[]
+function _printConcentrationsCSV(file::IOStream, simulation::Simulation{T}) where {T}
+    concentrations = getConcentrations(simulation.grid)
 
     for row in eachrow(concentrations)
         formatted_row = [Printf.@sprintf("%.6g", x) for x in row]  # Format each element like is done in the C++ version using Eigen3
         println(file, join(formatted_row, " "))
     end
-    println(file)  # Add extra newlines for separation
+    println(file)
     println(file)
 end
 
-function printConcentrations(simulation::Simulation{T,approach,solver}) where {T,approach,solver}
-    println(simulation.grid.concentrations[])
+function _printConcentrations(simulation::Simulation{T}) where {T}
+    println(getConcentrations(simulation.grid))
 end
 
-function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver}
+function run(simulation::Simulation{T}) where {T}
     file = nothing
     try
         if simulation.csvOutput > CSV_OUTPUT_OFF
-            file = createCSVfile(simulation)
+            file = _createCSVfile(simulation)
         end
 
         function simulationStepCallback()
             if simulation.consoleOutput >= CONSOLE_OUTPUT_VERBOSE
-                printConcentrations(simulation)
+                _printConcentrations(simulation)
             end
             if simulation.csvOutput >= CSV_OUTPUT_VERBOSE
-                printConcentrationsCSV(simulation, file)
+                _printConcentrationsCSV(file, simulation)
             end
         end
 
@@ -106,11 +110,11 @@ function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver
         end
 
         if simulation.consoleOutput >= CONSOLE_OUTPUT_ON
-            printConcentrations(simulation)
+            _printConcentrations(simulation)
         end
 
         if simulation.csvOutput >= CSV_OUTPUT_ON
-            printConcentrationsCSV(simulation, file)
+            _printConcentrationsCSV(file, simulation)
         end
 
     finally
@@ -120,18 +124,18 @@ function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver
     end
 end
 
-function setTimestep(simulation::Simulation{T,approach,solver}, timestep::T) where {T,approach,solver}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.iterations, timestep, simulation.consoleOutput, simulation.csvOutput)
+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
 
-function setIterations(simulation::Simulation{T,approach,solver}, iterations::Int) where {T,approach,solver}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, iterations, simulation.timestep, simulation.consoleOutput, simulation.csvOutput)
+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
 
-function setOutputConsole(simulation::Simulation{T,approach,solver}, consoleOutput::CONSOLE_OUTPUT) where {T,approach,solver}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.iterations, simulation.timestep, consoleOutput, simulation.csvOutput)
+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
 
-function setOutputCSV(simulation::Simulation{T,approach,solver}, csvOutput::CSV_OUTPUT) where {T,approach,solver}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.iterations, simulation.timestep, simulation.consoleOutput, csvOutput)
+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