diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index d0d7ff26..652dff84 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -5,6 +5,7 @@ using BenchmarkTools
 using ForwardDiff
 using IntervalArithmetic
 using IntervalRootFinding
+using StaticArrays
 
 import Random
 
@@ -43,6 +44,40 @@ for method in (Newton, Krawczyk)
 end
 
 
+S = SUITE["Linear equations"] = BenchmarkGroup()
+
+sizes = (2, 5, 10)
+
+for n in sizes
+    s = S["n = $n"] = BenchmarkGroup()
+    A = Interval.(randn(n, n))
+    b = Interval.(randn(n))
+    A = SMatrix{n, n}(A)
+    b = SVector{n}(b)
+
+    s["Gauss seidel"] = @benchmarkable gauss_seidel_interval($A, $b)
+    s["Gauss seidel contractor"] = @benchmarkable gauss_seidel_contractor($A, $b)
+    s["Gauss elimination"] = @benchmarkable gauss_elimination_interval($A, $b)
+end
+
+
+S = SUITE["Dietmar-Ratz"] = BenchmarkGroup()
+X = Interval(0.75, 1.75)
+
+for (k, dr) in enumerate(dr_functions)
+    s = S["Dietmar-Ratz $k"] = BenchmarkGroup()
+
+    if k != 8  # dr8 is excluded as it has too many roots
+        for method in (Newton, Krawczyk)
+            s[string(method)] = @benchmarkable roots($dr, $X, $method, $tol)
+        end
+    end
+
+    s["Automatic differentiation"] = @benchmarkable ForwardDiff.derivative($dr, $X)
+    s["Slope expansion"] = @benchmarkable slope($dr, $X, $(mid(X)))
+end
+
+
 S = SUITE["Linear equations"] = BenchmarkGroup()
 
 sizes = (2, 5, 10)
diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl
index d70ed783..3b80cae1 100644
--- a/src/IntervalRootFinding.jl
+++ b/src/IntervalRootFinding.jl
@@ -8,7 +8,7 @@ using IntervalArithmetic
 using ForwardDiff
 using StaticArrays
 
-using LinearAlgebra: I, Diagonal
+using LinearAlgebra: I, diag
 
 
 import Base: ⊆, show, big, \
@@ -38,19 +38,20 @@ const D = derivative
 
 const where_bisect = IntervalArithmetic.where_bisect ## 127//256
 
+include("complex.jl")
+include("slopes.jl")
+include("linear_eq.jl")
+include("quadratic.jl")
 
 include("root_object.jl")
 
-include("newton.jl")
-include("krawczyk.jl")
-
-include("complex.jl")
 include("contractors.jl")
 include("roots.jl")
+
+# Old stuff
+include("newton.jl")
+include("krawczyk.jl")
 include("newton1d.jl")
-include("quadratic.jl")
-include("linear_eq.jl")
-include("slopes.jl")
 
 
 gradient(f) = X -> ForwardDiff.gradient(f, X[:])
diff --git a/src/linear_eq.jl b/src/linear_eq.jl
index 1dc43b8a..39c14a2c 100644
--- a/src/linear_eq.jl
+++ b/src/linear_eq.jl
@@ -2,12 +2,10 @@
 Preconditions the matrix A and b with the inverse of mid(A)
 """
 function preconditioner(A::AbstractMatrix, b::AbstractArray)
-
     Aᶜ = mid.(A)
     B = inv(Aᶜ)
 
     return B*A, B*b
-
 end
 
 function gauss_seidel_interval(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
@@ -48,119 +46,103 @@ function gauss_seidel_interval!(x::AbstractArray, A::AbstractMatrix, b::Abstract
     x
 end
 
-function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+function gauss_seidel_contractor(A::SMatrix{N, N, Interval{T}},
+        b::SVector{N, Interval{T}} ;
+        precondition=true, maxiter=100) where {N, T}
+
+    x = MVector{N, Interval{T}}([entireinterval(T) for _ in 1:N])
+    Am = convert(MMatrix{N, N, Interval{T}}, A)
+    bm = convert(MVector{N, Interval{T}}, b)
+
+    gauss_seidel_contractor!(x, Am, bm, precondition=precondition, maxiter=maxiter)
 
-    n = size(A, 1)
-    x = similar(b)
-    x .= -1e16..1e16
-    x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter)
     return x
 end
 
-function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+function gauss_seidel_contractor!(x::MVector{N, Interval{T}},
+        A::MMatrix{N, N, Interval{T}},
+        b::MVector{N, Interval{T}} ;
+        precondition=true, maxiter=100) where {N, T}
 
     precondition && ((A, b) = preconditioner(A, b))
 
-    n = size(A, 1)
-
-    diagA = Diagonal(A)
+    diagA = diag(A)
     extdiagA = copy(A)
-    for i in 1:n
-        if (typeof(b) <: SArray)
-            extdiagA = setindex(extdiagA, Interval(0), i, i)
-        else
-            extdiagA[i, i] = Interval(0)
-        end
+    for i in 1:N
+        extdiagA[i, i] = Interval(zero(T))
     end
-    inv_diagA = inv(diagA)
+    inv_diagA = 1 ./ diagA
 
     for iter in 1:maxiter
         x¹ = copy(x)
-        x = x .∩ (inv_diagA * (b - extdiagA * x))
-        if all(x .== x¹)
-            break
-        end
+        x = x .∩ (inv_diagA .* (b - extdiagA * x))
+        all(x .== x¹) && return
     end
-    x
 end
 
-function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true)
+"""
+    gauss_elimination_interval(A::SMatrix, b::SVector, [precondition=true])
+"""
+function gauss_elimination_interval(A::SMatrix{N, N, Interval{T}},
+        b::SVector{N, Interval{T}} ;
+        precondition=true) where {N, T}
+
+    x = MVector{N, Interval{T}}([entireinterval(T) for _ in 1:N])
+    Am = convert(MMatrix{N, N, Interval{T}}, A)
+    bm = convert(MVector{N, Interval{T}}, b)
 
-    x = similar(b)
-    x .= -1e16..1e16
-    x = gauss_elimination_interval!(x, A, b, precondition=precondition)
+    gauss_elimination_interval!(x, Am, bm, precondition=precondition)
 
     return x
 end
+
 """
-Solves the system of linear equations using Gaussian Elimination.
+    gauss_elimination_interval!(x::MVector, A::MMatrix, b::MVector, [precondition=true])
+
+Solves the system of linear equations `A*x = b` using Gaussian Elimination.
 Preconditioning is used when the `precondition` keyword argument is `true`.
 
+This is the mutable version of the algorithm, which mutate `x`, `A` and `b`.
+Final upper triangle of `A` and `b` correspond to the triangular system.
+
 REF: Luc Jaulin et al.,
 *Applied Interval Analysis*, pg. 72
 """
-function gauss_elimination_interval!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true)
+function gauss_elimination_interval!(x::MVector{N, Interval{T}},
+        A::MMatrix{N, N, Interval{T}},
+        b::MVector{N, Interval{T}} ;
+        precondition=true) where {N, T}
 
     if precondition
-        (A, b) = preconditioner(A, b)
-    else
-        A = copy(A)
-        b = copy(b)
+        A, b = preconditioner(A, b)
     end
 
-    n = size(A, 1)
-
-    p = similar(b)
-    p .= 0
-
-    for i in 1:(n-1)
-        if 0 ∈ A[i, i] # diagonal matrix is not invertible
-            p .= entireinterval(b[1])
-            return p .∩ x  # return x?
+    for i in 1:(N-1)
+        if 0 ∈ A[i, i] # Pivot is zero, implying singular matrix
+            x .= entireinterval(T)
+            return
         end
 
-        for j in (i+1):n
+        for j in (i+1):N
             α = A[j, i] / A[i, i]
             b[j] -= α * b[i]
 
-            for k in (i+1):n
+            for k in (i+1):N
                 A[j, k] -= α * A[i, k]
             end
         end
     end
 
-    for i in n:-1:1
+    for i in N:-1:1
+        temp = zero(T)
 
-        temp = zero(b[1])
-
-        for j in (i+1):n
-            temp += A[i, j] * p[j]
+        for j in (i+1):N
+            temp += A[i, j] * x[j]
         end
 
-        p[i] = (b[i] - temp) / A[i, i]
+        x[i] = (b[i] - temp) / A[i, i]
     end
-
-    return p .∩ x
-end
-
-function gauss_elimination_interval1(A::AbstractMatrix, b::AbstractArray; precondition=true)
-
-    n = size(A, 1)
-    x = fill(-1e16..1e16, n)
-
-    x = gauss_elimination_interval1!(x, A, b, precondition=precondition)
-
-    return x
-end
-"""
-Using `Base.\``
-"""
-function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true)
-
-    precondition && ((a, b) = preconditioner(a, b))
-
-    a \ b
 end
 
-\(A::StaticMatrix{Interval{T}}, b::StaticArray{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)
-\(A::Matrix{Interval{T}}, b::Array{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)
+\(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}} ; kwargs...) where {N, T} = gauss_elimination_interval(A, b, kwargs...)
+\(A::SMatrix{N, N, Interval{T}}, b::IntervalBox{N, T} ; kwargs...) where {N, T} = gauss_elimination_interval(A, b.v, kwargs...)
diff --git a/test/linear_eq.jl b/test/linear_eq.jl
index 8b720998..358fb307 100644
--- a/test/linear_eq.jl
+++ b/test/linear_eq.jl
@@ -19,13 +19,13 @@ end
 
 @testset "Linear Equations" begin
 
-    As = [[2..3 0..1; 1..2 2..3], ]
-    bs = [[0..120, 60..240], ]
-    xs = [[-120..90, -60..240], ]
+    As = Any[SMatrix{2}(2..3, 1..2, 0..1, 2..3), ]
+    bs = Any[SVector(0..120, 60..240), ]
+    xs = Any[SVector(-120..90, -60..240), ]
 
-    for i in 1:10
-        rand_A = rand_mat(i)[1]
-        rand_x = rand_vec(i)[1]
+    for i in 2:6
+        rand_A = rand_mat(i)[3]
+        rand_x = rand_vec(i)[3]
         rand_b = rand_A * rand_x
         push!(As, rand_A)
         push!(bs, rand_b)
@@ -35,7 +35,7 @@ end
     n = length(As)
 
 
-    for solver in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval, \)
+    for solver in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval)
         @testset "Solver $solver" begin
             #for precondition in (false, true)
 
diff --git a/test/test_smiley.jl b/test/test_smiley.jl
index 19c8138c..889753ba 100644
--- a/test/test_smiley.jl
+++ b/test/test_smiley.jl
@@ -12,26 +12,28 @@ function test_all_unique(xs)
 end
 
 const tol = 1e-6
-const method = Newton # NOTE: Bisection method performs badly in all examples
-
-@info("Testing method $(method)")
+# NOTE: Bisection method performs badly in all examples
 
 @testset "$(SmileyExample22.title)" begin
-    roots_found = roots(SmileyExample22.f, SmileyExample22.region, method, tol)
-    @test length(roots_found) == 8
-    test_all_unique(roots_found)
-    # no reference data for roots given
+    for method in (Newton, Krawczyk)
+        roots_found = roots(SmileyExample22.f, SmileyExample22.region, method, tol)
+        @test length(roots_found) == 8
+        test_all_unique(roots_found)
+        # no reference data for roots given
+    end
 end
 
 for example in (SmileyExample52, SmileyExample54) #, SmileyExample55)
     @testset "$(example.title)" begin
-        roots_found = roots(example.f, example.region, method, tol)
-        @test length(roots_found) == length(example.known_roots)
-        test_all_unique(roots_found)
-        for rf in roots_found
-            # check there is exactly one known root for each found root
-            @test sum(!isempty(rk ∩ rf.interval)
-                    for rk in example.known_roots) == 1
+        for method in (Newton, Krawczyk)
+            roots_found = roots(example.f, example.region, method, tol)
+            @test length(roots_found) == length(example.known_roots)
+            test_all_unique(roots_found)
+            for rf in roots_found
+                # check there is exactly one known root for each found root
+                @test sum(!isempty(rk ∩ rf.interval)
+                        for rk in example.known_roots) == 1
+            end
         end
     end
 end