diff --git a/julia/tests/compare_csv.py b/julia/tests/compare_csv.py
new file mode 100644
index 0000000..f4fc179
--- /dev/null
+++ b/julia/tests/compare_csv.py
@@ -0,0 +1,82 @@
+import argparse
+import csv
+import sys
+
+
+def read_matrices_from_csv(file_path):
+    with open(file_path, newline='') as csvfile:
+        reader = csv.reader(csvfile, delimiter=' ', skipinitialspace=True)
+        matrices, matrix, line_numbers = [], [], []
+        line_number = 0
+
+        for row in reader:
+            line_number += 1
+            if row:
+                matrix.append([float(item) for item in row])
+                line_numbers.append(line_number)
+            else:
+                if matrix:
+                    matrices.append((matrix, line_numbers))
+                    matrix, line_numbers = [], []
+
+        if matrix:  # Add the last matrix if there's no blank line at the end
+            matrices.append((matrix, line_numbers))
+
+        return matrices
+
+def compare_matrices(matrix_a_tuple, matrix_b_tuple, tolerance=0.01):
+    data_a, line_numbers_a = matrix_a_tuple
+    data_b, line_numbers_b = matrix_b_tuple
+    max_diff = 0
+
+    if len(data_a) != len(data_b):
+        return False, f"Matrix size mismatch: Matrix A has {len(data_a)} rows, Matrix B has {len(data_b)} rows.", None, None, max_diff
+
+    for row_index, (row_a, row_b, line_num_a, line_num_b) in enumerate(zip(data_a, data_b, line_numbers_a, line_numbers_b)):
+        if len(row_a) != len(row_b):
+            return False, f"Row size mismatch at line {line_num_a} (File A) and line {line_num_b} (File B).", line_num_a, line_num_b, max_diff
+
+        for col_index, (elem_a, elem_b) in enumerate(zip(row_a, row_b)):
+            if abs(elem_a - elem_b) > max_diff:
+                max_diff = abs(elem_a - elem_b)
+            if abs(elem_a - elem_b) > tolerance:
+                return False, f"Mismatch at line {line_num_a}, column {col_index+1} (File A) vs line {line_num_b}, column {col_index+1} (File B):\nA({elem_a}) vs B({elem_b}).", line_num_a, line_num_b, max_diff
+
+    return True, "Matrices are equal.", None, None, max_diff
+
+def compare_csv_files(file_path1, file_path2, tolerance=0.01):
+    matrices1 = read_matrices_from_csv(file_path1)
+    matrices2 = read_matrices_from_csv(file_path2)
+    max_difference = 0
+
+    if len(matrices1) != len(matrices2):
+        return False, "Number of matrices in files do not match.", None, None, max_difference
+
+    for index, (matrix_a_tuple, matrix_b_tuple) in enumerate(zip(matrices1, matrices2)):
+        equal, message, line_num_a, line_num_b, max_diff = compare_matrices(matrix_a_tuple, matrix_b_tuple, tolerance)
+        if max_diff > max_difference:
+            max_difference = max_diff
+        if not equal:
+            return False, f"In Matrix pair {index}: {message}", line_num_a, line_num_b, max_difference
+
+    return True, "All matrices are equal.", None, None, max_difference
+
+
+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('--silent', action='store_true', help='Run in silent mode without printing details')
+    args = parser.parse_args()
+
+    are_equal, message, line_num_a, line_num_b, max_difference = compare_csv_files(args.file_a, args.file_b, args.tolerance)
+
+    if not args.silent:
+        print(message)
+
+    if not are_equal:
+        sys.exit(1)
+
+if __name__ == "__main__":
+    main()
diff --git a/julia/tests/cpp_bench/BTCS_100_100_1000.cpp b/julia/tests/cpp_bench/BTCS_100_100_1000.cpp
new file mode 100644
index 0000000..632c5a0
--- /dev/null
+++ b/julia/tests/cpp_bench/BTCS_100_100_1000.cpp
@@ -0,0 +1,39 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[])
+{
+    // **** GRID ****
+    int rows = 100;
+    int cols = 100;
+    Grid64 grid(rows, cols);
+
+    MatrixXd concentrations = MatrixXd::Constant(rows, cols, 0);
+    concentrations(10, 10) = 2000;
+    concentrations(90, 90) = 2000;
+    grid.setConcentrations(concentrations);
+
+    MatrixXd alphax = MatrixXd::Constant(rows, cols, 1);
+    MatrixXd alphay = MatrixXd::Constant(rows, cols, 1);
+    grid.setAlpha(alphax, alphay);
+
+    // **** BOUNDARY ****
+    Boundary bc = Boundary(grid);
+    bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
+    bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
+    bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
+    bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 2);
+
+    // **** SIMULATION ****
+    Simulation simulation = Simulation(grid, bc);
+    simulation.setTimestep(0.05);
+    simulation.setIterations(1000);
+    simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+    simulation.setOutputConsole(CONSOLE_OUTPUT_OFF);
+
+    // **** RUN SIMULATION ****
+    simulation.run();
+}
diff --git a/julia/tests/cpp_bench/BTCS_1_20_100.cpp b/julia/tests/cpp_bench/BTCS_1_20_100.cpp
new file mode 100644
index 0000000..651c5bd
--- /dev/null
+++ b/julia/tests/cpp_bench/BTCS_1_20_100.cpp
@@ -0,0 +1,34 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[])
+{
+  // **** GRID ****
+  int cells = 20;
+  Grid64 grid(cells);
+
+  MatrixXd concentrations = MatrixXd::Constant(1, cells, 0);
+  concentrations(0, 0) = 2000;
+  grid.setConcentrations(concentrations);
+
+  MatrixXd alpha = MatrixXd::Constant(1, cells, 1);
+  grid.setAlpha(alpha);
+
+  // **** BOUNDARY ****
+  Boundary bc = Boundary(grid);
+  bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
+  bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+
+  // **** SIMULATION ****
+  Simulation simulation = Simulation(grid, bc);
+  simulation.setTimestep(0.1);
+  simulation.setIterations(100);
+  simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+  simulation.setOutputConsole(CONSOLE_OUTPUT_OFF);
+
+  // **** RUN SIMULATION ****
+  simulation.run();
+}
diff --git a/julia/tests/cpp_bench/BTCS_1_45_750.cpp b/julia/tests/cpp_bench/BTCS_1_45_750.cpp
new file mode 100644
index 0000000..b98b1f8
--- /dev/null
+++ b/julia/tests/cpp_bench/BTCS_1_45_750.cpp
@@ -0,0 +1,36 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[])
+{
+    // **** GRID ****
+    int cells = 45;
+    Grid64 grid(cells);
+
+    MatrixXd concentrations = MatrixXd::Constant(1, cells, 10);
+    concentrations(0, 5) = 2000;
+    grid.setConcentrations(concentrations);
+
+    MatrixXd alpha = MatrixXd::Constant(1, cells, 1);
+    alpha.block(0, 0, 1, 15) = MatrixXd::Constant(1, 15, 0.5);
+    alpha.block(0, 30, 1, 15) = MatrixXd::Constant(1, 15, 1.5);
+    grid.setAlpha(alpha);
+
+    // **** BOUNDARY ****
+    Boundary bc = Boundary(grid);
+    bc.setBoundarySideClosed(BC_SIDE_LEFT);
+    bc.setBoundarySideClosed(BC_SIDE_RIGHT);
+
+    // **** SIMULATION ****
+    Simulation simulation = Simulation(grid, bc);
+    simulation.setTimestep(1.23);
+    simulation.setIterations(750);
+    simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+    simulation.setOutputConsole(CONSOLE_OUTPUT_OFF);
+
+    // **** RUN SIMULATION ****
+    simulation.run();
+}
diff --git a/julia/tests/cpp_bench/BTCS_20_20_500.cpp b/julia/tests/cpp_bench/BTCS_20_20_500.cpp
new file mode 100644
index 0000000..a509880
--- /dev/null
+++ b/julia/tests/cpp_bench/BTCS_20_20_500.cpp
@@ -0,0 +1,38 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[])
+{
+  // **** GRID ****
+  int rows = 20;
+  int cols = 20;
+  Grid64 grid(rows, cols);
+
+  MatrixXd concentrations = MatrixXd::Constant(rows, cols, 0);
+  concentrations(10, 10) = 2000;
+  grid.setConcentrations(concentrations);
+
+  MatrixXd alphax = MatrixXd::Constant(rows, cols, 1);
+  MatrixXd alphay = MatrixXd::Constant(rows, cols, 1);
+  grid.setAlpha(alphax, alphay);
+
+  // **** BOUNDARY ****
+  Boundary bc = Boundary(grid);
+  bc.setBoundarySideClosed(BC_SIDE_LEFT);
+  bc.setBoundarySideClosed(BC_SIDE_RIGHT);
+  bc.setBoundarySideClosed(BC_SIDE_TOP);
+  bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
+
+  // **** SIMULATION ****
+  Simulation simulation = Simulation(grid, bc);
+  simulation.setTimestep(0.1);
+  simulation.setIterations(500);
+  simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+  simulation.setOutputConsole(CONSOLE_OUTPUT_OFF);
+
+  // **** RUN SIMULATION ****
+  simulation.run();
+}
diff --git a/julia/tests/cpp_bench/BTCS_450_670_100.cpp b/julia/tests/cpp_bench/BTCS_450_670_100.cpp
new file mode 100644
index 0000000..963bd21
--- /dev/null
+++ b/julia/tests/cpp_bench/BTCS_450_670_100.cpp
@@ -0,0 +1,46 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[])
+{
+    // **** GRID ****
+    int rows = 450;
+    int cols = 670;
+    Grid64 grid(rows, cols);
+
+    MatrixXd concentrations = MatrixXd::Constant(rows, cols, 0);
+    concentrations(10, 10) = 1500;
+    concentrations(440, 660) = 750;
+    concentrations(440, 10) = 750;
+    concentrations(10, 660) = 750;
+    concentrations(220, 335) = 1500;
+    grid.setConcentrations(concentrations);
+
+    MatrixXd alphax = MatrixXd::Constant(rows, cols, 1);
+    MatrixXd alphay = MatrixXd::Constant(rows, cols, 1);
+    alphax.block(0, 0, 100, cols) = MatrixXd::Constant(100, cols, 0.5);
+    alphax.block(100, 0, 100, cols) = MatrixXd::Constant(100, cols, 0.8);
+    alphay.block(0, 0, rows, 200) = MatrixXd::Constant(rows, 200, 0.6);
+    alphay.block(0, 200, rows, 200) = MatrixXd::Constant(rows, 200, 0.9);
+    grid.setAlpha(alphax, alphay);
+
+    // **** BOUNDARY ****
+    Boundary bc = Boundary(grid);
+    bc.setBoundarySideClosed(BC_SIDE_LEFT);
+    bc.setBoundarySideClosed(BC_SIDE_RIGHT);
+    bc.setBoundarySideClosed(BC_SIDE_TOP);
+    bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
+
+    // **** SIMULATION ****
+    Simulation simulation = Simulation(grid, bc);
+    simulation.setTimestep(0.2);
+    simulation.setIterations(100);
+    simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+    simulation.setOutputConsole(CONSOLE_OUTPUT_OFF);
+
+    // **** RUN SIMULATION ****
+    simulation.run();
+}
diff --git a/julia/tests/julia_bench/BTCS_100_100_1000.jl b/julia/tests/julia_bench/BTCS_100_100_1000.jl
new file mode 100644
index 0000000..ed818e8
--- /dev/null
+++ b/julia/tests/julia_bench/BTCS_100_100_1000.jl
@@ -0,0 +1,36 @@
+include("../../tug/Simulation.jl")
+
+function main()
+    # **** GRID ****
+    rows::Int = 100
+    cols::Int = 100
+    grid::Grid = Grid{Float64}(rows, cols)
+
+    concentrations = fill(0.0, rows, cols)
+    concentrations[11, 11] = 2000
+    concentrations[91, 91] = 2000
+    setConcentrations!(grid, concentrations)
+
+    alphaX = fill(1.0, rows, cols)
+    alphaY = fill(1.0, rows, cols)
+    setAlpha!(grid, alphaX, alphaY)
+
+    # **** BOUNDARY ****
+    bc::Boundary = Boundary(grid)
+    setBoundarySideConstant!(bc, LEFT, 1.0)
+    setBoundarySideConstant!(bc, RIGHT, 1.0)
+    setBoundarySideConstant!(bc, TOP, 0.0)
+    setBoundarySideConstant!(bc, BOTTOM, 2.0)
+
+    # **** SIMULATION ****
+    simulation::Simulation = Simulation(grid, bc)
+    simulation = setTimestep(simulation, 0.05)
+    simulation = setIterations(simulation, 1000)
+    simulation = setOutputConsole(simulation, CONSOLE_OUTPUT_OFF)
+    simulation = setOutputCSV(simulation, CSV_OUPUT_VERBOSE)
+
+    # **** RUN SIMULATION ****
+    run(simulation)
+end
+
+main()
diff --git a/julia/tests/julia_bench/BTCS_1_20_100.jl b/julia/tests/julia_bench/BTCS_1_20_100.jl
new file mode 100644
index 0000000..2267425
--- /dev/null
+++ b/julia/tests/julia_bench/BTCS_1_20_100.jl
@@ -0,0 +1,31 @@
+include("../../tug/Simulation.jl")
+
+function main()
+    # **** GRID ****
+    cells::Int = 20
+    grid::Grid = Grid{Float64}(cells)
+
+    concentrations = fill(0.0, 1, cells)
+    concentrations[1] = 2000
+    setConcentrations!(grid, concentrations)
+
+    alpha = fill(1.0, 1, cells)
+    setAlpha!(grid, alpha)
+
+    # **** BOUNDARY ****
+    bc::Boundary = Boundary(grid)
+    setBoundarySideConstant!(bc, LEFT, 0.0)
+    setBoundarySideConstant!(bc, RIGHT, 0.0)
+
+    # **** SIMULATION ****
+    simulation::Simulation = Simulation(grid, bc)
+    simulation = setTimestep(simulation, 0.1)
+    simulation = setIterations(simulation, 100)
+    simulation = setOutputConsole(simulation, CONSOLE_OUTPUT_OFF)
+    simulation = setOutputCSV(simulation, CSV_OUPUT_VERBOSE)
+
+    # **** RUN SIMULATION ****
+    run(simulation)
+end
+
+main()
diff --git a/julia/tests/julia_bench/BTCS_1_45_750.jl b/julia/tests/julia_bench/BTCS_1_45_750.jl
new file mode 100644
index 0000000..5f38221
--- /dev/null
+++ b/julia/tests/julia_bench/BTCS_1_45_750.jl
@@ -0,0 +1,33 @@
+include("../../tug/Simulation.jl")
+
+function main()
+    # **** GRID ****
+    cells::Int = 45
+    grid::Grid = Grid{Float64}(cells)
+
+    concentrations = fill(10.0, 1, cells)
+    concentrations[6] = 2000
+    setConcentrations!(grid, concentrations)
+
+    alpha = fill(1.0, 1, cells)
+    alpha[1:15] .= 0.5
+    alpha[31:45] .= 1.5
+    setAlpha!(grid, alpha)
+
+    # **** BOUNDARY ****
+    bc::Boundary = Boundary(grid)
+    setBoundarySideClosed!(bc, LEFT)
+    setBoundarySideClosed!(bc, RIGHT)
+
+    # **** SIMULATION ****
+    simulation::Simulation = Simulation(grid, bc)
+    simulation = setTimestep(simulation, 1.23)
+    simulation = setIterations(simulation, 750)
+    simulation = setOutputConsole(simulation, CONSOLE_OUTPUT_OFF)
+    simulation = setOutputCSV(simulation, CSV_OUPUT_VERBOSE)
+
+    # **** RUN SIMULATION ****
+    run(simulation)
+end
+
+main()
diff --git a/julia/tests/julia_bench/BTCS_20_20_500.jl b/julia/tests/julia_bench/BTCS_20_20_500.jl
new file mode 100644
index 0000000..5ec0e18
--- /dev/null
+++ b/julia/tests/julia_bench/BTCS_20_20_500.jl
@@ -0,0 +1,35 @@
+include("../../tug/Simulation.jl")
+
+function main()
+    # **** GRID ****
+    rows::Int = 20
+    cols::Int = 20
+    grid::Grid = Grid{Float64}(rows, cols)
+
+    concentrations = fill(0.0, rows, cols)
+    concentrations[11, 11] = 2000
+    setConcentrations!(grid, concentrations)
+
+    alphaX = fill(1.0, rows, cols)
+    alphaY = fill(1.0, rows, cols)
+    setAlpha!(grid, alphaX, alphaY)
+
+    # **** BOUNDARY ****
+    bc::Boundary = Boundary(grid)
+    setBoundarySideClosed!(bc, LEFT)
+    setBoundarySideClosed!(bc, RIGHT)
+    setBoundarySideClosed!(bc, TOP)
+    setBoundarySideClosed!(bc, BOTTOM)
+
+    # **** SIMULATION ****
+    simulation::Simulation = Simulation(grid, bc)
+    simulation = setTimestep(simulation, 0.1)
+    simulation = setIterations(simulation, 500)
+    simulation = setOutputConsole(simulation, CONSOLE_OUTPUT_OFF)
+    simulation = setOutputCSV(simulation, CSV_OUPUT_VERBOSE)
+
+    # **** RUN SIMULATION ****
+    run(simulation)
+end
+
+main()
diff --git a/julia/tests/julia_bench/BTCS_450_670_100.jl b/julia/tests/julia_bench/BTCS_450_670_100.jl
new file mode 100644
index 0000000..cb2dfba
--- /dev/null
+++ b/julia/tests/julia_bench/BTCS_450_670_100.jl
@@ -0,0 +1,43 @@
+include("../../tug/Simulation.jl")
+
+function main()
+    # **** GRID ****
+    rows::Int = 450
+    cols::Int = 670
+    grid::Grid = Grid{Float64}(rows, cols)
+
+    concentrations = fill(0.0, rows, cols)
+    concentrations[11, 11] = 1500
+    concentrations[441, 661] = 750
+    concentrations[441, 11] = 750
+    concentrations[11, 661] = 750
+    concentrations[221, 336] = 1500
+    setConcentrations!(grid, concentrations)
+
+    alphaX = fill(1.0, rows, cols)
+    alphaY = fill(1.0, rows, cols)
+    alphaX[1:100, :] .= 0.5
+    alphaX[101:200, :] .= 0.8
+    alphaY[:, 1:200] .= 0.6
+    alphaY[:, 201:400] .= 0.9
+    setAlpha!(grid, alphaX, alphaY)
+
+    # **** BOUNDARY ****
+    bc::Boundary = Boundary(grid)
+    setBoundarySideClosed!(bc, LEFT)
+    setBoundarySideClosed!(bc, RIGHT)
+    setBoundarySideClosed!(bc, TOP)
+    setBoundarySideClosed!(bc, BOTTOM)
+
+    # **** SIMULATION ****
+    simulation::Simulation = Simulation(grid, bc)
+    simulation = setTimestep(simulation, 0.2)
+    simulation = setIterations(simulation, 100)
+    simulation = setOutputConsole(simulation, CONSOLE_OUTPUT_OFF)
+    simulation = setOutputCSV(simulation, CSV_OUPUT_VERBOSE)
+
+    # **** RUN SIMULATION ****
+    run(simulation)
+end
+
+main()
diff --git a/julia/tests/test.py b/julia/tests/test.py
new file mode 100644
index 0000000..0642ae3
--- /dev/null
+++ b/julia/tests/test.py
@@ -0,0 +1,131 @@
+import os
+import subprocess
+import argparse
+from compare_csv import compare_csv_files
+import time
+
+# ANSI color codes
+RED = '\033[0;31m'
+GREEN = '\033[0;32m'
+NC = '\033[0m'  # No Color
+
+def remove_non_empty_dir(path):
+    if os.path.exists(path):
+        for root, dirs, files in os.walk(path, topdown=False):
+            for name in files:
+                os.remove(os.path.join(root, name))
+            for name in dirs:
+                os.rmdir(os.path.join(root, name))
+        os.rmdir(path)
+
+def get_max_name_length(directory):
+    max_length = 0
+    for file in os.listdir(directory):
+        if file.endswith('.csv'):
+            name_length = len(os.path.splitext(file)[0])
+            if name_length > max_length:
+                max_length = name_length
+    return max_length
+
+def format_difference(diff):
+    return '{:.5f}'.format(diff).rjust(8) if diff != 0 else '0'.rjust(8)
+
+def run_benchmark(command, runs):
+    total_time = 0
+    for _ in range(runs):
+        start_time = time.time()
+        subprocess.run(command)
+        total_time += time.time() - start_time
+    return total_time / runs
+
+def main(tolerance, runs, silent):
+    BENCHMARK_DIR = "./cpp_bench"
+    JULIA_DIR = "./julia_bench"
+    COMPILER = "g++"
+    CFLAGS = ["-O3", "-fopenmp", "-I", "/usr/local/include", "-I", "../../include/", "-I", "/usr/include/eigen3"]
+    BIN_DIR = "./cpp_bin_temp"
+    OUTPUT_DIR = "./csv_temp"
+
+    # Clean up and create directories
+    for dir_path in [BIN_DIR, OUTPUT_DIR]:
+        remove_non_empty_dir(dir_path)
+        os.makedirs(dir_path)
+
+    for file in os.listdir('.'):
+        if file.endswith('.csv'):
+            os.remove(file)
+
+    # Compile and run C++ benchmarks
+    if not silent: print("-----  Running C++ Benchmarks  -----")
+    cpp_times = {}
+    for benchmark in os.listdir(BENCHMARK_DIR):
+        if benchmark.endswith(".cpp"):
+            name = os.path.splitext(benchmark)[0]
+            if not silent: print(f"Compiling {name}...", end="", flush=True)
+            subprocess.run([COMPILER, *CFLAGS, "-o", f"{BIN_DIR}/{name}", f"{BENCHMARK_DIR}/{benchmark}"])
+            if not silent: print(" Running...", end="", flush=True)
+            cpp_times[name] = run_benchmark([f"./{BIN_DIR}/{name}"], runs)
+            if not silent: print(" Done.", flush=True)
+
+    # Move CSV files to output directory
+    for file in os.listdir('.'):
+        if file.endswith('.csv'):
+            os.rename(file, f"{OUTPUT_DIR}/{file}")
+
+    max_name_length = get_max_name_length(OUTPUT_DIR)
+
+    # Run Julia benchmarks and compare
+    if not silent: print("\n----- Running Julia Benchmarks -----")
+    results_dict = {}
+    pass_all = True
+    julia_times = {}
+    for csv_file in sorted(os.listdir(OUTPUT_DIR), key=lambda x: os.path.splitext(x)[0]):
+        name = os.path.splitext(csv_file)[0]
+        padded_name = name.ljust(max_name_length)
+        if os.path.exists(f"{JULIA_DIR}/{name}.jl"):
+            if not silent: print(f"Running {name}...", end="", flush=True)
+            julia_times[name] = run_benchmark(["julia", f"{JULIA_DIR}/{name}.jl"], runs)
+
+            if os.path.exists(f"./{name}.csv"):
+                are_equal, _, _, _, max_diff = compare_csv_files(f"./{name}.csv", f"{OUTPUT_DIR}/{csv_file}", tolerance)
+                formatted_diff = format_difference(max_diff)
+                cpp_time = '{:.4f}s'.format(cpp_times[name]).rjust(8)
+                julia_time = '{:.4f}s'.format(julia_times[name]).rjust(8)
+                result = f"{padded_name}: {'Success' if are_equal else 'Failure'} (Max Diff: {formatted_diff}, C++: {cpp_time}, Julia: {julia_time})"
+                result_color = GREEN if are_equal else RED
+                results_dict[name] = f"{result_color}{result}{NC}"
+                if not are_equal:
+                    pass_all = False
+
+            else:
+                results_dict[name] = f"{RED}{padded_name}: No Julia output{NC}"
+                pass_all = False
+            if not silent: print(" Done.", flush=True)
+            
+
+    # Clean up
+    remove_non_empty_dir(BIN_DIR)
+    remove_non_empty_dir(OUTPUT_DIR)
+    
+    for file in os.listdir('.'):
+        if file.endswith('.csv'):
+            os.remove(file)
+
+    # Print results
+    if not silent: print("\n-----     Benchmark Results    -----")
+    print(f"Parameters: Tolerance = {tolerance}, Runs = {runs}")
+    for name in sorted(results_dict):
+        print(results_dict[name])
+    if pass_all:
+        print(f"\n{GREEN}All benchmarks and comparisons passed.{NC}")
+    else:
+        print(f"\n{RED}Some benchmarks or comparisons failed.{NC}")
+        exit(1)
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description='Benchmark and Compare Script')
+    parser.add_argument('--tolerance', type=float, default=0.005, help='Tolerance for CSV comparison')
+    parser.add_argument('--runs', type=int, default=1, help='Number of benchmark runs')
+    parser.add_argument('--silent', action='store_true', help='Run in silent mode without printing details')
+    args = parser.parse_args()
+    main(args.tolerance, args.runs, args.silent)
diff --git a/julia/tug/Boundary.jl b/julia/tug/Boundary.jl
new file mode 100644
index 0000000..455f7e8
--- /dev/null
+++ b/julia/tug/Boundary.jl
@@ -0,0 +1,102 @@
+# Boundary.jl
+# Julia implementation of Boundary types, translating from C++'s Boundary.hpp
+
+@enum TYPE CLOSED CONSTANT
+@enum SIDE LEFT = 1 RIGHT = 2 TOP = 3 BOTTOM = 4
+
+
+# BoundaryElement class
+struct BoundaryElement{T}
+    type::TYPE
+    value::T
+
+    # Generic constructor for closed case
+    BoundaryElement{T}() where {T} = new{T}(CLOSED, convert(T, -1))
+
+    # Constructor for constant case
+    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
+
+function getValue(be::BoundaryElement{T})::T where {T}
+    be.value
+end
+
+function setValue!(be::BoundaryElement{T}, value::T) where {T}
+    if value < 0
+        throw(ArgumentError("No negative concentration allowed."))
+    end
+    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
+    cols::UInt32
+    rows::UInt32
+    boundaries::Vector{Vector{BoundaryElement{T}}}
+
+    # Constructor
+    function Boundary(grid::Grid{T}) where {T}
+        dim = grid.dim
+        cols = grid.cols
+        rows = grid.rows
+        boundaries = Vector{Vector{BoundaryElement{T}}}(undef, 4)
+
+        if dim == 1
+            boundaries[Int(LEFT)] = [BoundaryElement{T}()]
+            boundaries[Int(RIGHT)] = [BoundaryElement{T}()]
+        elseif dim == 2
+            boundaries[Int(LEFT)] = [BoundaryElement{T}() for _ in 1:rows]
+            boundaries[Int(RIGHT)] = [BoundaryElement{T}() for _ in 1:rows]
+            boundaries[Int(TOP)] = [BoundaryElement{T}() for _ in 1:cols]
+            boundaries[Int(BOTTOM)] = [BoundaryElement{T}() for _ in 1:cols]
+        else
+            throw(ArgumentError("Only 1- and 2-dimensional grids are defined!"))
+        end
+
+        new{T}(dim, cols, rows, boundaries)
+    end
+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."))
+    end
+
+    is_vertical = side in (LEFT, RIGHT)
+    n = is_vertical ? boundary.rows : boundary.cols
+
+    boundary.boundaries[Int(side)] = [BoundaryElement{T}() for _ in 1:n]
+end
+
+function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::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
+
+    is_vertical = side in (LEFT, RIGHT)
+    n = is_vertical ? boundary.rows : boundary.cols
+
+    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
new file mode 100644
index 0000000..1046b77
--- /dev/null
+++ b/julia/tug/Core/BTCS.jl
@@ -0,0 +1,244 @@
+# BTCS.jl
+# Implementation of heterogenous BTCS (backward time-centered space)
+# solution of diffusion equation in 1D and 2D space using the
+# alternating-direction implicit (ADI) method.
+
+using LinearAlgebra
+using SparseArrays
+
+include("../Boundary.jl")
+include("../Grid.jl")
+
+# Helper functions and types
+function calcAlphaIntercell(alpha1::T, alpha2::T, useHarmonic::Bool=true) where {T}
+    if useHarmonic
+        return 2 / ((1 / alpha1) + (1 / alpha2))
+    else
+        return 0.5 * (alpha1 + alpha2)
+    end
+end
+
+# calculates coefficient for boundary in constant case
+function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where {T}
+    centerCoeff = 1 + sx * (calcAlphaIntercell(alpha_center, alpha_side) + 2 * alpha_center)
+    sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side)
+    return (centerCoeff, sideCoeff)
+end
+
+# calculates coefficient for boundary in closed case
+function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T}
+    centerCoeff = 1 + sx * calcAlphaIntercell(alpha_center, alpha_side)
+    sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side)
+    return (centerCoeff, sideCoeff)
+end
+
+# creates coefficient matrix for next time step from alphas in x-direction
+function createCoeffMatrix(alpha::Matrix{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, numCols::Int, rowIndex::Int, sx::T) where {T}
+    numCols = max(numCols, 2)
+    cm = spzeros(T, numCols, numCols)
+
+    # left column
+    if getType(bcLeft[rowIndex]) == CONSTANT
+        centerCoeffTop, rightCoeffTop = calcBoundaryCoeffConstant(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
+        cm[1, 1] = centerCoeffTop
+        cm[1, 2] = rightCoeffTop
+    elseif getType(bcLeft[rowIndex]) == CLOSED
+        centerCoeffTop, rightCoeffTop = calcBoundaryCoeffClosed(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
+        cm[1, 1] = centerCoeffTop
+        cm[1, 2] = rightCoeffTop
+    else
+        error("Undefined Boundary Condition Type somewhere on Left or Top!")
+    end
+
+    # inner columns
+    for i in 2:(numCols-1)
+        cm[i, i-1] = -sx * calcAlphaIntercell(alpha[rowIndex, i-1], alpha[rowIndex, i])
+        cm[i, i] = 1 + sx * (calcAlphaIntercell(alpha[rowIndex, i], alpha[rowIndex, i+1]) + calcAlphaIntercell(alpha[rowIndex, i-1], alpha[rowIndex, i]))
+        cm[i, i+1] = -sx * calcAlphaIntercell(alpha[rowIndex, i], alpha[rowIndex, i+1])
+    end
+
+    # right column
+    if getType(bcRight[rowIndex]) == CONSTANT
+        centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeffConstant(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
+        cm[numCols, numCols-1] = leftCoeffBottom
+        cm[numCols, numCols] = centerCoeffBottom
+    elseif getType(bcRight[rowIndex]) == CLOSED
+        centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeffClosed(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
+        cm[numCols, numCols-1] = leftCoeffBottom
+        cm[numCols, numCols] = centerCoeffBottom
+    else
+        error("Undefined Boundary Condition Type somewhere on Right or Bottom!")
+    end
+
+    return cm
+end
+
+
+
+
+function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
+    sy * calcAlphaIntercell(alpha_center, alpha_neighbor) * conc_center +
+    (1 - sy * calcAlphaIntercell(alpha_center, alpha_neighbor)) * conc_center
+end
+
+
+function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
+    sy * calcAlphaIntercell(alpha_center, alpha_neighbor) * conc_center +
+    (1 - sy * (calcAlphaIntercell(alpha_center, alpha_center) + 2 * alpha_center)) * conc_center +
+    sy * alpha_center * conc_bc
+end
+
+
+function createSolutionVector(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}}, length::Int, rowIndex::Int, sx::T, sy::T) where {T}
+    numRows = size(concentrations, 1)
+    sv = Vector{T}(undef, length)
+
+    # Inner rows
+    if rowIndex > 1 && rowIndex < numRows
+        for i = 1:length
+            sv[i] = sy * calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i]) * concentrations[rowIndex+1, i] +
+                    (1 - sy * (calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i]) + calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]))) * concentrations[rowIndex, i] +
+                    sy * calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) * concentrations[rowIndex-1, i]
+        end
+    end
+
+    # First row
+    if rowIndex == 1
+        for i = 1:length
+            if getType(bcTop[i]) == CONSTANT
+                sv[i] = calcExplicitConcentrationsBoundaryConstant(concentrations[rowIndex, i], getValue(bcTop[i]), alphaY[rowIndex, i], alphaY[rowIndex+1, i], sy)
+            elseif getType(bcTop[i]) == CLOSED
+                sv[i] = calcExplicitConcentrationsBoundaryClosed(concentrations[rowIndex, i], alphaY[rowIndex, i], alphaY[rowIndex+1, i], sy)
+            else
+                error("Undefined Boundary Condition Type somewhere on Left or Top!")
+            end
+        end
+    end
+
+    # Last row
+    if rowIndex == numRows
+        for i = 1:length
+            if getType(bcBottom[i]) == CONSTANT
+                sv[i] = calcExplicitConcentrationsBoundaryConstant(concentrations[rowIndex, i], getValue(bcBottom[i]), alphaY[rowIndex, i], alphaY[rowIndex-1, i], sy)
+            elseif getType(bcBottom[i]) == CLOSED
+                sv[i] = calcExplicitConcentrationsBoundaryClosed(concentrations[rowIndex, i], alphaY[rowIndex, i], alphaY[rowIndex-1, i], sy)
+            else
+                error("Undefined Boundary Condition Type somewhere on Right or Bottom!")
+            end
+        end
+    end
+
+    # First column - additional fixed concentration change from perpendicular dimension in constant BC case
+    if getType(bcLeft[rowIndex]) == CONSTANT
+        sv[1] += 2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex])
+    end
+
+    # Last column - additional fixed concentration change from perpendicular dimension in constant BC case
+    if getType(bcRight[rowIndex]) == CONSTANT
+        sv[end] += 2 * sx * alphaX[rowIndex, end] * getValue(bcRight[rowIndex])
+    end
+
+    return sv
+end
+
+
+# solver for linear equation system; A corresponds to coefficient matrix, b to the solution vector
+function LinearAlgebraAlgorithm(A::SparseMatrixCSC{T}, b::Vector{T}) where {T}
+    return A \ b
+end
+
+# BTCS solution for 1D grid
+function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Function) where {T}
+    length = grid.cols
+    sx = timestep / (grid.deltaCol * grid.deltaCol)
+
+    b = Vector{T}(undef, length)
+
+    alpha = grid.alphaX[]
+
+    bcLeft = getBoundarySide(bc, LEFT)
+    bcRight = getBoundarySide(bc, RIGHT)
+
+    concentrations = grid.concentrations[]
+    rowIndex = 1
+    A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx)
+    for i in 1:length
+        b[i] = concentrations[1, i]
+    end
+    if getType(getBoundarySide(bc, LEFT)[1]) == CONSTANT
+        b[1] += 2 * sx * alpha[1, 1] * bcLeft[1].value
+    end
+    if getType(getBoundarySide(bc, RIGHT)[1]) == CONSTANT
+        b[length] += 2 * sx * alpha[1, length] * bcRight[1].value
+    end
+
+    concentrations_t1 = solverFunc(A, b)
+
+    for j in 1:length
+        concentrations[1, j] = concentrations_t1[j]
+    end
+
+    setConcentrations!(grid, concentrations)
+end
+
+# BTCS solution for 2D grid
+function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Function, numThreads::Int) where {T}
+    rowMax = grid.rows
+    colMax = grid.cols
+    sx = timestep / (2 * grid.deltaCol * grid.deltaCol)
+    sy = timestep / (2 * grid.deltaRow * grid.deltaRow)
+
+    A = spzeros(T, rowMax, rowMax)
+    concentrations_t1 = zeros(T, rowMax, colMax)
+    row_t1 = Vector{T}(undef, colMax)
+
+    alphaX = grid.alphaX[]
+    alphaY = grid.alphaY[]
+
+    bcLeft = getBoundarySide(bc, LEFT)
+    bcRight = getBoundarySide(bc, RIGHT)
+    bcTop = getBoundarySide(bc, TOP)
+    bcBottom = getBoundarySide(bc, BOTTOM)
+
+    concentrations = grid.concentrations[]
+
+    for i = 1:rowMax
+        A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx)
+        b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy)
+
+        row_t1 = solverFunc(A, b)
+
+        concentrations_t1[i, :] = row_t1
+    end
+
+    concentrations_t1 = copy(transpose(concentrations_t1))
+    concentrations = copy(transpose(concentrations))
+    alphaX = copy(transpose(alphaX))
+    alphaY = copy(transpose(alphaY))
+
+    for i = 1:colMax
+        # Swap alphas, boundary conditions and sx/sy for column-wise calculation
+        A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy)
+        b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, bcLeft, bcRight, rowMax, i, sy, sx)
+
+        row_t1 = solverFunc(A, b)
+
+        concentrations[i, :] = row_t1
+    end
+
+    concentrations = copy(transpose(concentrations))
+
+    setConcentrations!(grid, concentrations)
+end
+
+
+# Entry point for EigenLU solver; differentiate between 1D and 2D grid
+function BTCS_LU(grid::Grid{T}, bc::Boundary{T}, timestep::T, numThreads::Int=1) where {T}
+    if grid.dim == 1
+        BTCS_1D(grid, bc, timestep, LinearAlgebraAlgorithm)
+    elseif grid.dim == 2
+        BTCS_2D(grid, bc, timestep, LinearAlgebraAlgorithm, numThreads)
+    else
+        error("Error: Only 1- and 2-dimensional grids are defined!")
+    end
+end
diff --git a/julia/tug/Grid.jl b/julia/tug/Grid.jl
new file mode 100644
index 0000000..1c3d476
--- /dev/null
+++ b/julia/tug/Grid.jl
@@ -0,0 +1,59 @@
+using LinearAlgebra
+
+struct Grid{T}
+    cols::Int
+    rows::Int
+    dim::Int
+    domainCol::T
+    domainRow::T
+    deltaCol::T
+    deltaRow::T
+    concentrations::Ref{Matrix{T}}
+    alphaX::Ref{Matrix{T}}
+    alphaY::Ref{Matrix{T}}
+
+    # Constructor for 1D-Grid
+    function Grid{T}(length::Int) where {T}
+        if length <= 3
+            throw(ArgumentError("Given grid length too small. Must be greater than 3."))
+        end
+
+        new{T}(length, 1, 1, T(length), 0, T(1), 0, Ref(fill(T(0), 1, length)), Ref(fill(T(0), 1, length)), Ref(fill(T(0), 1, length)))
+    end
+
+    # Constructor for 2D-Grid
+    function Grid{T}(row::Int, col::Int) where {T}
+        if row <= 3 || col <= 3
+            throw(ArgumentError("Given grid dimensions too small. Must each be greater than 3."))
+        end
+
+        new{T}(col, row, 2, T(col), T(row), T(1), T(1), Ref(fill(T(0), row, col)), Ref(fill(T(0), row, col)), Ref(fill(T(0), row, col)))
+    end
+end
+
+function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where {T}
+    grid.concentrations[] = new_concentrations
+end
+
+function setAlpha!(grid::Grid{T}, alpha::Matrix{T}) where {T}
+    if grid.dim != 1
+        error("Grid is not one dimensional, you should probably use the 2D setter function!")
+    end
+    if size(alpha, 1) != 1 || size(alpha, 2) != grid.cols
+        error("Given matrix of alpha coefficients mismatch with Grid dimensions!")
+    end
+
+    grid.alphaX[] = alpha
+end
+
+function setAlpha!(grid::Grid{T}, alphaX::Matrix{T}, alphaY::Matrix{T}) where {T}
+    if grid.dim != 2
+        error("Grid is not two dimensional, you should probably use the 1D setter function!")
+    end
+    if size(alphaX) != (grid.rows, grid.cols) || size(alphaY) != (grid.rows, grid.cols)
+        error("Given matrices of alpha coefficients mismatch with Grid dimensions!")
+    end
+
+    grid.alphaX[] = alphaX
+    grid.alphaY[] = alphaY
+end
diff --git a/julia/tug/Simulation.jl b/julia/tug/Simulation.jl
new file mode 100644
index 0000000..70fdc9a
--- /dev/null
+++ b/julia/tug/Simulation.jl
@@ -0,0 +1,140 @@
+using Printf
+
+include("Grid.jl")
+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_OUPUT CSV_OUPUT_OFF CSV_OUPUT_ON CSV_OUPUT_VERBOSE CSV_OUPUT_XTREME
+
+# Create the Simulation class
+struct Simulation{T,approach,solver}
+    grid::Grid{T}
+    bc::Boundary{T}
+    approach::APPROACH
+    solver::SOLVER
+
+    innerIterations::Int
+    iterations::Int
+    numThreads::Int
+    timestep::T
+
+    consoleOutput::CONSOLE_OUTPUT
+    csvOutput::CSV_OUPUT
+
+    function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS,
+        solver::SOLVER=EIGEN_LU_SOLVER, innerIterations::Int=1, iterations::Int=1,
+        numThreads::Int=1, timestep::T=0.1,
+        consoleOutput::CONSOLE_OUTPUT=CONSOLE_OUTPUT_OFF, csvOutput::CSV_OUPUT=CSV_OUPUT_OFF) where {T}
+        new{T,APPROACH,SOLVER}(grid, bc, approach, solver, innerIterations, iterations, numThreads, timestep, consoleOutput, csvOutput)
+    end
+
+end
+
+
+function createCSVfile(simulation::Simulation{T,approach,solver}) 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
+    numIterations = simulation.iterations
+    filename = string(approachString, "_", row, "_", col, "_", numIterations, ".csv")
+
+    while isfile(filename)
+        appendIdent += 1
+        filename = string(approachString, "_", row, "_", col, "_", numIterations, "-", appendIdent, ".csv")
+    end
+
+    open(filename, "w") do file
+        # Write boundary conditions if required
+        if simulation.csvOutput == CSV_OUPUT_XTREME
+            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)
+            end
+
+            write(file, "\n\n")
+        end
+    end
+
+    return filename
+end
+
+function writeBoundarySideValues(file, 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}, filename::String) where {T,approach,solver}
+    concentrations = simulation.grid.concentrations[]
+
+    open(filename, "a") do file  # Open file in append mode
+        for row in eachrow(concentrations)
+            println(file, join(row, " "))
+        end
+        println(file)  # Add extra newlines for separation
+        println(file)
+    end
+end
+
+function printConcentrations(simulation::Simulation{T,approach,solver}) where {T,approach,solver}
+    println(simulation.grid.concentrations[])
+end
+
+function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver}
+    filename::String = ""
+    if simulation.csvOutput > CSV_OUPUT_OFF
+        filename = createCSVfile(simulation)
+    end
+
+    if simulation.approach == BTCS
+        if simulation.solver == EIGEN_LU_SOLVER
+            for i in 1:(simulation.iterations*simulation.innerIterations)
+                if simulation.consoleOutput == CONSOLE_OUTPUT_VERBOSE
+                    printConcentrations(simulation)
+                end
+
+                if simulation.csvOutput >= CSV_OUPUT_VERBOSE
+                    printConcentrationsCSV(simulation, filename)
+                end
+
+                BTCS_LU(simulation.grid, simulation.bc, simulation.timestep, simulation.numThreads)
+            end
+        else
+            error("Undefined solver!")
+        end
+    else
+        error("Undefined approach!")
+    end
+
+    if simulation.consoleOutput == CONSOLE_OUTPUT_ON || simulation.consoleOutput == CONSOLE_OUTPUT_VERBOSE
+        printConcentrations(simulation)
+    end
+
+    if simulation.csvOutput == CSV_OUPUT_ON || simulation.csvOutput == CSV_OUPUT_VERBOSE || simulation.csvOutput == CSV_OUPUT_XTREME
+        printConcentrationsCSV(simulation, filename)
+    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.innerIterations, simulation.iterations, simulation.numThreads, 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, simulation.innerIterations, iterations, simulation.numThreads, simulation.timestep, simulation.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.innerIterations, simulation.iterations, simulation.numThreads, simulation.timestep, consoleOutput, simulation.csvOutput)
+end
+
+function setOutputCSV(simulation::Simulation{T,approach,solver}, csvOutput::CSV_OUPUT) where {T,approach,solver}
+    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.innerIterations, simulation.iterations, simulation.numThreads, simulation.timestep, simulation.consoleOutput, csvOutput)
+end