diff --git a/Dionysos.jl b/Dionysos.jl
new file mode 160000
index 000000000..de849a593
--- /dev/null
+++ b/Dionysos.jl
@@ -0,0 +1 @@
+Subproject commit de849a593d1be30840e9be25521e563d1d7e82bd
diff --git a/test/runtests.jl b/test/runtests.jl
index beadad195..01a581287 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,6 +6,8 @@ include("./utils/rectangle.jl")
 include("./utils/ellipsoid.jl")
 include("./utils/scalar_functions.jl")
 include("./utils/test_lazy_set_operations.jl")
+include("./utils/optim/newton_method.jl")
+include("./utils/optim/bisection.jl")
 
 include("./domain/test_griddomain.jl")
 include("./domain/test_general_domain.jl")
diff --git a/test/utils/optim/bisection.jl b/test/utils/optim/bisection.jl
new file mode 100644
index 000000000..54b04e6eb
--- /dev/null
+++ b/test/utils/optim/bisection.jl
@@ -0,0 +1,29 @@
+using Test
+using Dionysos
+bisection = Dionysos.Utils.bisection
+
+# Define the test function, its derivative, and second derivative
+@testset "Bissection" begin
+    # Define a test function
+    test_function(x) = x^2 - 4
+    ϵ = 1e-6
+
+    # Call the bisection function
+    result = bisection(test_function; interval = [0, 4], δ = ϵ, verbose = false)
+    @test isapprox(result[1], -4.0; atol = ϵ) && 0 ≤ result[2] ≤ 4
+
+    result = bisection(test_function; interval = [0, 4], δ = ϵ, verbose = true)
+    @test isapprox(result[1], -4.0; atol = ϵ) && 0 ≤ result[2] ≤ 4
+
+    result = bisection(test_function; interval = [1, 3], δ = ϵ, verbose = true)
+    @test isapprox(result[1], -3.0; atol = ϵ) && 0 ≤ result[2] ≤ 4
+
+    result = bisection(
+        test_function;
+        interval = [-5, 4],
+        δ = ϵ,
+        verbose = false,
+        stopIfNegative = true,
+    )
+    @test (-5 ≤ result[2] ≤ 4) && abs(result[1] + 4.0) > ϵ && result[1] < 0
+end
diff --git a/test/utils/optim/newton_method.jl b/test/utils/optim/newton_method.jl
new file mode 100644
index 000000000..a6efa531d
--- /dev/null
+++ b/test/utils/optim/newton_method.jl
@@ -0,0 +1,41 @@
+using Test
+using Dionysos
+newton_method = Dionysos.Utils.newton_method
+
+# Define the test function, its derivative, and second derivative
+@testset "EllipsoidBasics" begin
+    f(x) = x^2 - 4
+    df(x) = 2x
+    ddf(x) = 2
+
+    # Test parameters
+    interval = [1, 3]
+    x0 = 1.5
+    ϵ = 1e-6
+
+    # Call the newton_method function with the test parameters
+    result = newton_method(
+        f,
+        df,
+        ddf;
+        interval = interval,
+        x0 = x0,
+        ϵ = ϵ,
+        verbose = false,
+        stopIfNegative = false,
+    )
+    # Check the result
+    @test isapprox(result[1], -3.0; atol = ϵ) && interval[1] ≤ result[2] ≤ interval[2]
+
+    result = newton_method(
+        f,
+        df,
+        ddf;
+        interval = [0, 4],
+        x0 = x0,
+        ϵ = ϵ,
+        verbose = false,
+        stopIfNegative = false,
+    )
+    @test isapprox(result[1], -4.0; atol = ϵ) && 0 ≤ result[2] ≤ 4
+end