diff --git a/julia/tests/DistributedTest.jl b/julia/tests/DistributedTest.jl
new file mode 100644
index 0000000..39f9841
--- /dev/null
+++ b/julia/tests/DistributedTest.jl
@@ -0,0 +1,175 @@
+using BenchmarkTools
+using CSV
+using DataFrames
+using Distributed
+using Statistics
+
+include("../tug/TUG.jl")
+using .TUG
+
+# 1. Environment Setup
+function setup_environment(num_procs::Int)
+    if num_procs > 0
+        added_procs = addprocs(num_procs)
+
+        # Use remotecall to include TUG on each new worker
+        for proc in added_procs
+            remotecall_wait(include, proc, "../tug/TUG.jl")
+            remotecall_wait(eval, proc, :(using .TUG))
+        end
+    end
+end
+
+# 2. Test Case Definitions
+function testBTCS100()::Tuple{Grid,Boundary}
+    rows::Int = 1024
+    cols::Int = 1000
+
+    alphaX = fill(1.25, rows, cols)
+    alphaY = fill(1.1, rows, cols)
+    alphaX[1:100, :] .= 0.5
+    alphaX[101:200, :] .= 0.8
+    alphaY[:, 1:200] .= 0.6
+    alphaY[:, 201:400] .= 0.9
+    grid::Grid = Grid{Float64}(rows, cols, alphaX, alphaY)
+
+    concentrations = fill(0.5, rows, cols)
+    concentrations[11, 11] = 15000
+    concentrations[1015, 991] = 7500
+    concentrations[11, 991] = 7500
+    concentrations[1015, 11] = 7500
+    setConcentrations!(grid, concentrations)
+
+    bc::Boundary = Boundary(grid)
+    setBoundarySideClosed!(bc, LEFT)
+    setBoundarySideClosed!(bc, RIGHT)
+    setBoundarySideClosed!(bc, TOP)
+    setBoundarySideClosed!(bc, BOTTOM)
+
+    return grid, bc
+end
+
+function testBTCS200()::Tuple{Grid,Boundary}
+    rows::Int = 2027
+    cols::Int = 1999
+
+    alphaX = [sin(i / 100) * cos(j / 100) for i in 1:rows, j in 1:cols]
+    alphaY = [cos(i / 100) * sin(j / 100) for i in 1:rows, j in 1:cols]
+
+    grid::Grid = Grid{Float64}(rows, cols, alphaX, alphaY)
+
+    concentrations = [i * j / 1e2 for i in 1:rows, j in 1:cols]
+    concentrations[11, 11] = 15000
+    concentrations[2021, 1995] = 7500
+    concentrations[11, 1995] = 7500
+    concentrations[2021, 11] = 7500
+    setConcentrations!(grid, concentrations)
+
+    bc::Boundary = Boundary(grid)
+    setBoundarySideClosed!(bc, LEFT)
+    setBoundarySideConstant!(bc, RIGHT, 1.5)
+    setBoundarySideClosed!(bc, TOP)
+    setBoundarySideConstant!(bc, BOTTOM, 0.75)
+
+    return grid, bc
+end
+
+function testFTCS500()::Tuple{Grid,Boundary}
+    rows::Int = 2000
+    cols::Int = 2000
+
+    alphaX = [sin(i / 100) * cos(j / 100) + 1 for i in 1:rows, j in 1:cols]
+    alphaY = [cos(i / 100) * sin(j / 100) + 1 for i in 1:rows, j in 1:cols]
+
+    grid::Grid = Grid{Float64}(rows, cols, alphaX, alphaY)
+
+    concentrations = [(i * j) / 1e6 for i in 1:rows, j in 1:cols]
+    concentrations[1001, 1001] = 2000
+    setConcentrations!(grid, concentrations)
+
+    bc::Boundary = Boundary(grid)
+    setBoundarySideClosed!(bc, LEFT)
+    setBoundarySideClosed!(bc, RIGHT)
+    setBoundarySideClosed!(bc, TOP)
+    setBoundarySideConstant!(bc, BOTTOM, 0.75)
+
+    return grid, bc
+end
+
+# 3. Simulation Runners
+function run_static_simulation(grid::Grid, bc::Boundary, method, steps::Int, dt::Float64)
+    simulation = Simulation(grid, bc, method, steps, dt, CONSOLE_OUTPUT_OFF, CSV_OUTPUT_OFF)
+    TUG.run(simulation)
+end
+
+function run_dynamic_simulation(grid::Grid, bc::Boundary, method, steps::Int, dt::Float64, num_procs::Int)
+    setup_environment(num_procs)
+
+    simulation = DynamicSimulation(grid, bc, method, dt)
+    createGrid(simulation)
+    for _ in 1:steps
+        next(simulation)
+    end
+
+    if num_procs > 0
+        rmprocs(workers())
+    end
+end
+
+# 4. Performance Metrics
+function measure_performance(test_case_function, method, steps::Int, dt::Float64, num_procs::Int, dynamic::Bool=false)
+    grid, bc = test_case_function()
+
+    if dynamic
+        simulation_run = @benchmark run_dynamic_simulation($grid, $bc, $method, $steps, $dt, $num_procs)
+    else
+        simulation_run = @benchmark run_static_simulation($grid, $bc, $method, $steps, $dt)
+    end
+
+    time_measurement = mean(simulation_run).time / 1e9
+    memory_measurement = mean(simulation_run).memory / 1e6
+
+    return (time=time_measurement, memory=memory_measurement)
+end
+
+# 5. Benchmarking and Comparison
+function benchmark_test_case(test_case_function, method, steps::Int, dt::Float64, is_distributed::Bool, num_procs::Int)
+    performance_metrics = measure_performance(test_case_function, method, steps, dt, num_procs, is_distributed)
+    return performance_metrics
+end
+
+# 6. Reporting
+function report_results(results)
+    for result in results
+        println("Test Case: $(result.test_case)$(result.is_distributed ? ", Procs: $(result.num_procs)" :  ""), Time: $(result.time) s, Memory: $(result.memory) MB")
+    end
+end
+
+# 7. Cleanup
+function cleanup()
+    for filename in readdir()
+        if occursin(".csv", filename)
+            rm(filename, force=true)
+        end
+    end
+end
+
+# Main Function
+function main()
+    test_cases = [(testBTCS100, BTCS, 100, 0.01), (testBTCS200, BTCS, 200, 0.005), (testFTCS500, FTCS, 500, 0.005)]
+    configurations = [(false, 0), (true, 0), (true, 1), (true, 2)]  # Non-distributed, distributed with 1 and 2 additional procs
+
+    results = []
+
+    for (test_case, method, steps, dt) in test_cases
+        for (is_distributed, num_procs) in configurations
+            performance_metrics = benchmark_test_case(test_case, method, steps, dt, is_distributed, num_procs)
+            push!(results, (test_case=test_case, method=method, is_distributed=is_distributed, num_procs=num_procs, time=performance_metrics.time, memory=performance_metrics.memory))
+        end
+    end
+
+    cleanup()
+    report_results(results)
+end
+
+main()
diff --git a/julia/tests/compare_csv.py b/julia/tests/compare_csv.py
index f4fc179..4184d78 100644
--- a/julia/tests/compare_csv.py
+++ b/julia/tests/compare_csv.py
@@ -44,7 +44,7 @@ def compare_matrices(matrix_a_tuple, matrix_b_tuple, tolerance=0.01):
 
     return True, "Matrices are equal.", None, None, max_diff
 
-def compare_csv_files(file_path1, file_path2, tolerance=0.01):
+def compare_csv_files(file_path1, file_path2, tolerance=0):
     matrices1 = read_matrices_from_csv(file_path1)
     matrices2 = read_matrices_from_csv(file_path2)
     max_difference = 0
@@ -66,7 +66,7 @@ def main():
     parser = argparse.ArgumentParser(description='Compare two CSV files containing matrices.')
     parser.add_argument('file_a', help='Path to File A')
     parser.add_argument('file_b', help='Path to File B')
-    parser.add_argument('--tolerance', type=float, default=0.01, help='Tolerance for comparison (default: 0.01)')
+    parser.add_argument('--tolerance', type=float, default=0, help='Tolerance for comparison (default: 0)')
     parser.add_argument('--silent', action='store_true', help='Run in silent mode without printing details')
     args = parser.parse_args()
 
diff --git a/julia/tug/AbstractSimulation.jl b/julia/tug/AbstractSimulation.jl
new file mode 100644
index 0000000..7e3eb7d
--- /dev/null
+++ b/julia/tug/AbstractSimulation.jl
@@ -0,0 +1,67 @@
+using Printf
+
+include("Boundary.jl")
+include("Grid.jl")
+
+@enum APPROACH BTCS FTCS
+@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
+
+abstract type AbstractSimulation{T} end
+
+function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing)::IOStream where {T}
+    grid = grid === nothing ? simulation.grid : grid
+    approachString = string(simulation.approach)
+    rows = getRows(grid)
+    cols = getCols(grid)
+    numIterations = simulation.iterations
+
+    filename = string(approachString, "_", rows, "_", cols, "_", numIterations, ".csv")
+
+    appendIdent = 0
+    while isfile(filename)
+        appendIdent += 1
+        filename = string(approachString, "_", rows, "_", cols, "_", numIterations, "-", appendIdent, ".csv")
+    end
+
+    # Write boundary conditions if required
+    if simulation.csvOutput >= CSV_OUTPUT_XTREME
+        open(filename, "w") do file
+            writeBoundarySideValuesToCSV(file, simulation.bc, LEFT)
+            writeBoundarySideValuesToCSV(file, simulation.bc, RIGHT)
+
+            if getDim(grid) == 2
+                writeBoundarySideValuesToCSV(file, simulation.bc, TOP)
+                writeBoundarySideValuesToCSV(file, simulation.bc, BOTTOM)
+            end
+
+            write(file, "\n\n")
+        end
+    end
+
+    file = open(filename, "a")
+    return file
+end
+
+function writeBoundarySideValuesToCSV(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 writeConcentrationsToCLI(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing) where {T}
+    grid = grid === nothing ? simulation.grid : grid
+    println(getConcentrations(grid))
+end
+
+function writeConcentrationsToCSV(file::IOStream, simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing) where {T}
+    grid = grid === nothing ? simulation.grid : grid
+    concentrations = getConcentrations(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)
+    println(file)
+end
diff --git a/julia/tug/Boundary.jl b/julia/tug/Boundary.jl
index 88810e3..834ced1 100644
--- a/julia/tug/Boundary.jl
+++ b/julia/tug/Boundary.jl
@@ -3,10 +3,11 @@
 # condition at the edges of the diffusion grid. 
 # Translated from C++'s Boundary.hpp.
 
+include("Grid.jl")
+
 @enum TYPE CLOSED CONSTANT
 @enum SIDE LEFT = 1 RIGHT = 2 TOP = 3 BOTTOM = 4
 
-
 # BoundaryElement class
 struct BoundaryElement{T}
     type::TYPE
diff --git a/julia/tug/Core/BTCS.jl b/julia/tug/Core/BTCS.jl
index 16b9058..fa47a4b 100644
--- a/julia/tug/Core/BTCS.jl
+++ b/julia/tug/Core/BTCS.jl
@@ -4,9 +4,9 @@
 # alternating-direction implicit (ADI) method.
 # Translated from C++'s BTCS.hpp.
 
+using Base.Threads
 using LinearAlgebra
 using SparseArrays
-using Base.Threads
 
 include("../Boundary.jl")
 include("../Grid.jl")
diff --git a/julia/tug/Core/FTCS.jl b/julia/tug/Core/FTCS.jl
index 4afd953..c4d28b9 100644
--- a/julia/tug/Core/FTCS.jl
+++ b/julia/tug/Core/FTCS.jl
@@ -3,6 +3,10 @@
 # solution of diffusion equation in 1D and 2D space.
 # Translated from C++'s FTCS.hpp. 
 
+using Base.Threads
+using LinearAlgebra
+using SparseArrays
+
 include("../Boundary.jl")
 include("../Grid.jl")
 include("Utils.jl")
diff --git a/julia/tug/DynamicSimulation.jl b/julia/tug/DynamicSimulation.jl
new file mode 100644
index 0000000..a7a6336
--- /dev/null
+++ b/julia/tug/DynamicSimulation.jl
@@ -0,0 +1,69 @@
+using Distributed
+
+include("AbstractSimulation.jl")
+include("Boundary.jl")
+include("Core/BTCS.jl")
+include("Core/FTCS.jl")
+include("Grid.jl")
+
+struct DynamicSimulation{T} <: AbstractSimulation{T}
+    grid::Grid{T}
+    grids::Vector{Grid{T}}
+    bc::Boundary{T}
+
+    approach::APPROACH
+    iterations::Int
+    timestep::T
+
+    # Constructor
+    function DynamicSimulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH, timestep::T) where {T}
+        timestep, iterations = adjustTimestep(grid, approach, timestep, 1)
+        new{T}(grid, Vector{Grid{T}}(), bc, approach, iterations, timestep)
+    end
+end
+
+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}
+    pmap(grid -> runSimulationForGrid(simulation, grid), simulation.grids)
+end
+
+function printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int) where {T}
+    writeConcentrationsToCLI(simulation, simulation.grids[gridIndex])
+end
+
+function printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int) 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}
+    if simulation.approach == BTCS
+        runBTCS(grid, simulation.bc, simulation.timestep, simulation.iterations, () -> nothing)
+    elseif simulation.approach == FTCS
+        runFTCS(grid, simulation.bc, simulation.timestep, simulation.iterations, () -> nothing)
+    else
+        error("Undefined approach!")
+    end
+end
+
+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.grids[gridIndex], concentrations)
+end
+
+function setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T}) where {T}
+    setAlphaX!(simulation.grids[gridIndex], alphaX)
+end
+
+function setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T}) where {T}
+    setAlphaY!(simulation.grids[gridIndex], alphaY)
+end
diff --git a/julia/tug/Grid.jl b/julia/tug/Grid.jl
index 884bf9c..8390511 100644
--- a/julia/tug/Grid.jl
+++ b/julia/tug/Grid.jl
@@ -5,7 +5,7 @@
 
 using LinearAlgebra
 
-# Grid class
+# Grid class 
 struct Grid{T}
     cols::Int
     rows::Int
@@ -15,10 +15,10 @@ struct Grid{T}
     deltaCol::T
     deltaRow::T
     concentrations::Ref{Matrix{T}}
-    alphaX::Matrix{T}
-    alphaY::Union{Matrix{T},Nothing}
-    alphaX_t::Union{Matrix{T},Nothing}
-    alphaY_t::Union{Matrix{T},Nothing}
+    alphaX::Ref{Matrix{T}}
+    alphaY::Union{Ref{Matrix{T}},Nothing}
+    alphaX_t::Union{Ref{Matrix{T}},Nothing}
+    alphaY_t::Union{Ref{Matrix{T}},Nothing}
 
     # Constructor for 1D-Grid
     function Grid{T}(length::Int, alpha::Matrix{T}) where {T}
@@ -47,22 +47,30 @@ struct Grid{T}
 
         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}
+        new{T}(cols, rows, dim, domainCol, domainRow, deltaCol, deltaRow, concentrations, alphaX, alphaY, alphaX_t, alphaY_t)
+    end
+end
+
+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
 
 function getAlphaX(grid::Grid{T})::Matrix{T} where {T}
-    grid.alphaX
+    grid.alphaX[]
 end
 
 function getAlphaY(grid::Grid{T})::Matrix{T} where {T}
-    grid.alphaY
+    grid.alphaY[]
 end
 
 function getAlphaX_t(grid::Grid{T})::Matrix{T} where {T}
-    grid.alphaX_t
+    grid.alphaX_t[]
 end
 
 function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
-    grid.alphaY_t
+    grid.alphaY_t[]
 end
 
 function getCols(grid::Grid{T})::Int where {T}
@@ -89,6 +97,16 @@ function getRows(grid::Grid{T})::Int where {T}
     grid.rows
 end
 
+function setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T}) where {T}
+    grid.alphaX[] = new_alphaX
+    grid.alphaX_t[] = new_alphaX'
+end
+
+function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T}) where {T}
+    grid.alphaY[] = new_alphaY
+    grid.alphaY_t[] = new_alphaY'
+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 264ec25..1239a47 100644
--- a/julia/tug/Simulation.jl
+++ b/julia/tug/Simulation.jl
@@ -4,19 +4,14 @@
 # options. Simulation object also holds a predefined Grid and Boundary object.
 # Translated from C++'s Simulation.hpp.
 
-using Printf
-
-include("Grid.jl")
+include("AbstractSimulation.jl")
 include("Boundary.jl")
 include("Core/BTCS.jl")
 include("Core/FTCS.jl")
-
-@enum APPROACH BTCS FTCS
-@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
+include("Grid.jl")
 
 # Simulation class
-struct Simulation{T}
+struct Simulation{T} <: AbstractSimulation{T}
     grid::Grid{T}
     bc::Boundary{T}
     approach::APPROACH
@@ -34,7 +29,7 @@ struct Simulation{T}
     end
 end
 
-function _adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iterations::Int)::Tuple{T,Int} where {T}
+function adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iterations::Int)::Tuple{T,Int} where {T}
     if approach == FTCS
         if getDim(grid) == 1
             deltaSquare = getDeltaCol(grid)
@@ -64,91 +59,37 @@ function _adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iterati
     return timestep, iterations
 end
 
-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, "_", rows, "_", cols, "_", numIterations, ".csv")
-
-    appendIdent = 0
-    while isfile(filename)
-        appendIdent += 1
-        filename = string(approachString, "_", rows, "_", cols, "_", numIterations, "-", appendIdent, ".csv")
-    end
-
-    # Write boundary conditions if required
-    if simulation.csvOutput >= CSV_OUTPUT_XTREME
-        open(filename, "w") do file
-            _writeBoundarySideValues(file, simulation.bc, LEFT)
-            _writeBoundarySideValues(file, simulation.bc, RIGHT)
-
-            if getDim(simulation.grid) == 2
-                _writeBoundarySideValues(file, simulation.bc, TOP)
-                _writeBoundarySideValues(file, simulation.bc, BOTTOM)
-            end
-
-            write(file, "\n\n")
-        end
-    end
-
-    file = open(filename, "a")
-    return file
-end
-
-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(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)
-    println(file)
-end
-
-function _printConcentrations(simulation::Simulation{T}) where {T}
-    println(getConcentrations(simulation.grid))
-end
-
 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)
+                writeConcentrationsToCLI(simulation)
             end
             if simulation.csvOutput >= CSV_OUTPUT_VERBOSE
-                _printConcentrationsCSV(file, simulation)
+                writeConcentrationsToCSV(file, simulation)
             end
         end
 
         if simulation.approach == BTCS
             runBTCS(simulation.grid, simulation.bc, simulation.timestep, simulation.iterations, simulationStepCallback)
         elseif simulation.approach == FTCS
-            timestep, iterations = _adjustTimestep(simulation.grid, simulation.approach, simulation.timestep, simulation.iterations)
+            timestep, iterations = adjustTimestep(simulation.grid, simulation.approach, simulation.timestep, simulation.iterations)
             runFTCS(simulation.grid, simulation.bc, timestep, iterations, simulationStepCallback)
         else
             error("Undefined approach!")
         end
 
         if simulation.consoleOutput >= CONSOLE_OUTPUT_ON
-            _printConcentrations(simulation)
+            writeConcentrationsToCLI(simulation)
         end
 
         if simulation.csvOutput >= CSV_OUTPUT_ON
-            _printConcentrationsCSV(file, simulation)
+            writeConcentrationsToCSV(file, simulation)
         end
 
     finally
diff --git a/julia/tug/TUG.jl b/julia/tug/TUG.jl
new file mode 100644
index 0000000..93b8388
--- /dev/null
+++ b/julia/tug/TUG.jl
@@ -0,0 +1,30 @@
+module TUG
+
+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
+
+include("Simulation.jl")
+
+export Simulation
+export run, setTimestep, setIterations, setOutputConsole, setOutputCSV
+
+include("DynamicSimulation.jl")
+
+export DynamicSimulation
+export createGrid, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!, next, printConcentrationsCSV, printConcentrations
+
+include("Boundary.jl")
+
+export Boundary, BoundaryElement, TYPE, SIDE, BC_TYPE_CLOSED, BC_TYPE_CONSTANT, LEFT, RIGHT, TOP, BOTTOM
+export getType, getValue, setValue!, setBoundarySideClosed!, setBoundarySideConstant!, getBoundarySide, getBoundarySideClosed, getBoundarySideConstant
+
+include("Grid.jl")
+
+export Grid
+export getAlphaX, getAlphaY, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!, getDomainCol, getDomainRow, getDeltaCol, getDeltaRow, getDim
+
+include("Core/BTCS.jl")
+include("Core/FTCS.jl")
+
+end