diff --git a/examples/linear_eq.jl b/examples/linear_eq.jl
new file mode 100644
index 00000000..918bcb69
--- /dev/null
+++ b/examples/linear_eq.jl
@@ -0,0 +1,32 @@
+# Examples from Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis
+
+using IntervalArithmetic, IntervalRootFinding, StaticArrays
+
+A = [4..5 -1..1 1.5..2.5; -0.5..0.5 -7.. -5 1..2; -1.5.. -0.5 -0.7.. -0.5 2..3]
+sA = SMatrix{3}{3}(A)
+mA = MMatrix{3}{3}(A)
+
+b = [3..4, 0..2, 3..4]
+sb = SVector{3}(b)
+mb = MVector{3}(b)
+
+p = fill(-1e16..1e16, 3)
+
+rts = gauss_seidel_interval!(p, A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval!(p, sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval!(p, mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+
+rts = gauss_seidel_interval(A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval(sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval(mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+
+rts = gauss_seidel_contractor!(p, A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor!(p, sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor!(p, mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+
+rts = gauss_seidel_contractor(A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor(sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor(mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+
+rts = gauss_elimination_interval!(p, A, b, precondition=true) # Gaussian Elimination; precondition=true by default
+rts = gauss_elimination_interval(A, b, precondition=true) # Gaussian Elimination; precondition=true by default
diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl
new file mode 100644
index 00000000..a418bfba
--- /dev/null
+++ b/perf/linear_eq.jl
@@ -0,0 +1,63 @@
+using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays
+
+function randVec(n::Int)
+    a = randn(n)
+    A = Interval.(a)
+    mA = MVector{n}(A)
+    sA = SVector{n}(A)
+    return A, mA, sA
+end
+
+function randMat(n::Int)
+    a = randn(n, n)
+    A = Interval.(a)
+    mA = MMatrix{n, n}(A)
+    sA = SMatrix{n, n}(A)
+    return A, mA, sA
+end
+
+function benchmark(max=10)
+    df = DataFrame()
+    df[:Method] = ["Array", "MArray", "SArray", "Contractor", "ContractorMArray", "ContractorSArray"]
+    for n in 1:max
+        A, mA, sA = randMat(n)
+        b, mb, sb = randVec(n)
+        t1 = @belapsed gauss_seidel_interval($A, $b)
+        t2 = @belapsed gauss_seidel_interval($mA, $mb)
+        t3 = @belapsed gauss_seidel_interval($sA, $sb)
+        t4 = @belapsed gauss_seidel_contractor($A, $b)
+        t5 = @belapsed gauss_seidel_contractor($mA, $mb)
+        t6 = @belapsed gauss_seidel_contractor($sA, $sb)
+        df[Symbol("$n")] = [t1, t2, t3, t4, t5, t6]
+    end
+    a = []
+    for i in 1:max
+        push!(a, Symbol("$i"))
+    end
+    df1 = stack(df, a)
+    dfnew = unstack(df1, :variable, :Method, :value)
+    dfnew = rename(dfnew, :variable => :n)
+    println(dfnew)
+    dfnew
+end
+
+function benchmark_elimination(max=10)
+    df = DataFrame()
+    df[:Method] = ["Gauss Elimination", "Base.\\"]
+    for n in 1:max
+        A, mA, sA = randMat(n)
+        b, mb, sb = randVec(n)
+        t1 = @belapsed gauss_elimination_interval($A, $b)
+        t2 = @belapsed gauss_elimination_interval1($A, $b)
+        df[Symbol("$n")] = [t1, t2]
+    end
+    a = []
+    for i in 1:max
+        push!(a, Symbol("$i"))
+    end
+    df1 = stack(df, a)
+    dfnew = unstack(df1, :variable, :Method, :value)
+    dfnew = rename(dfnew, :variable => :n)
+    println(dfnew)
+    dfnew
+end
diff --git a/perf/linear_eq_results.txt b/perf/linear_eq_results.txt
new file mode 100644
index 00000000..ff08cefe
--- /dev/null
+++ b/perf/linear_eq_results.txt
@@ -0,0 +1,149 @@
+For n = 1
+  11.382 μs (420 allocations: 16.83 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-0.858728, -0.858727]]
+  10.980 μs (421 allocations: 16.77 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-0.858728, -0.858727]]
+
+For n = 2
+  27.322 μs (822 allocations: 33.84 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[0.946446, 0.946447], [0.40483, 0.404831]]
+  26.037 μs (821 allocations: 32.41 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[0.946446, 0.946447], [0.40483, 0.404831]]
+
+For n = 3
+  50.094 μs (1222 allocations: 50.20 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[1.21978, 1.21979], [-2.29245, -2.29244], [-5.4913, -5.49129]]
+  47.146 μs (1221 allocations: 48.06 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[1.21978, 1.21979], [-2.29245, -2.29244], [-5.4913, -5.49129]]
+
+For n = 4
+  78.244 μs (1626 allocations: 66.92 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-1.23305, -1.23304], [-0.267412, -0.267411], [-27.3277, -27.3276], [-8.3822, -8.38219]]
+  74.046 μs (1621 allocations: 63.70 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-1.23305, -1.23304], [-0.267412, -0.267411], [-27.3277, -27.3276], [-8.3822, -8.38219]]
+
+For n = 5
+  109.184 μs (2026 allocations: 83.63 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[0.670434, 0.670435], [0.876143, 0.876144], [0.0960248, 0.0960249], [-0.109509, -0.109508], [0.375848, 0.375849]]
+  112.329 μs (2032 allocations: 83.41 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[0.670434, 0.670435], [0.876143, 0.876144], [0.0960248, 0.0960249], [-0.109509, -0.109508], [0.375848, 0.375849]]
+
+For n = 6
+  148.001 μs (2426 allocations: 100.09 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-2.40426, -2.40425], [-2.04742, -2.04741], [-1.39128, -1.39127], [-1.28434, -1.28433], [3.76269, 3.7627], [-1.28934, -1.28933]]
+  148.599 μs (2432 allocations: 99.73 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-2.40426, -2.40425], [-2.04742, -2.04741], [-1.39128, -1.39127], [-1.28434, -1.28433], [3.76269, 3.7627], [-1.28934, -1.28933]]
+
+For n = 7
+  191.907 μs (2826 allocations: 116.92 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[0.706613, 0.706614], [2.97926, 2.97927], [-2.55176, -2.55175], [-1.51898, -1.51897], [-0.987214, -0.987213], [-2.07295, -2.07294], [2.29913, 2.29914]]
+  194.024 μs (2832 allocations: 116.34 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[0.706613, 0.706614], [2.97926, 2.97927], [-2.55176, -2.55175], [-1.51898, -1.51897], [-0.987214, -0.987213], [-2.07295, -2.07294], [2.29913, 2.29914]]
+
+For n = 8
+  246.387 μs (3226 allocations: 133.58 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-0.149055, -0.149054], [7.1198, 7.11981], [-8.30751, -8.3075], [-2.37305, -2.37304], [0.0150316, 0.0150317], [3.35756, 3.35757], [3.9637, 3.96371], [-4.7502, -4.75019]]
+  248.388 μs (3232 allocations: 132.72 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-0.149055, -0.149054], [7.1198, 7.11981], [-8.30751, -8.3075], [-2.37305, -2.37304], [0.0150316, 0.0150317], [3.35756, 3.35757],[3.9637, 3.96371], [-4.7502, -4.75019]]
+
+For n = 9
+  313.552 μs (3626 allocations: 150.41 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-0.319338, -0.319337], [1.05902, 1.05903], [-1.87228, -1.87227], [-0.878659, -0.878658], [1.2596, 1.25961], [-2.90817, -2.90816],[-2.77684, -2.77683], [-3.76999, -3.76998], [-0.355478, -0.355477]]
+  309.842 μs (3632 allocations: 149.23 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-0.319338, -0.319337], [1.05902, 1.05903], [-1.87228, -1.87227], [-0.878659, -0.878658], [1.2596, 1.25961], [-2.90817, -2.90816], [-2.77684, -2.77683], [-3.76999, -3.76998], [-0.355478, -0.355477]]
+
+For n = 10
+  383.613 μs (4026 allocations: 167.34 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[2.98301, 2.98302], [0.995662, 0.995663], [0.150873, 0.150874], [1.03224, 1.03225], [-0.887891, -0.88789], [0.547691, 0.547692], [-0.530325, -0.530324], [5.39277, 5.39278], [0.230642, 0.230643], [-1.87602, -1.87601]]
+  393.352 μs (4032 allocations: 165.84 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[2.98301, 2.98302], [0.995662, 0.995663], [0.150873, 0.150874], [1.03224, 1.03225], [-0.887891, -0.88789], [0.547691, 0.547692], [-0.530325, -0.530324], [5.39277, 5.39278], [0.230642, 0.230643], [-1.87602, -1.87601]]
+
+For n = 11
+  460.411 μs (4426 allocations: 184.31 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-1.49959, -1.49958], [0.280678, 0.280679], [-1.53757, -1.53756], [1.2517, 1.25171], [-0.605238, -0.605237], [0.825546, 0.825547],[0.755109, 0.75511], [0.18521, 0.185211], [-1.82804, -1.82803], [0.384323, 0.384324], [0.241048, 0.241049]]
+  458.773 μs (4432 allocations: 182.59 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-1.49959, -1.49958], [0.280678, 0.280679], [-1.53757, -1.53756], [1.2517, 1.25171], [-0.605238, -0.605237], [0.825546, 0.825547], [0.755109, 0.75511], [0.18521, 0.185211], [-1.82804, -1.82803], [0.384323, 0.384324], [0.241048, 0.241049]]
+
+For n = 12
+  550.373 μs (4826 allocations: 201.45 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-2.79797, -2.79796], [-1.73355, -1.73354], [-0.606508, -0.606507], [-2.44883, -2.44882], [1.71346, 1.71347], [-2.72265, -2.72264], [-1.55253, -1.55252], [0.289974, 0.289975], [0.252046, 0.252047], [0.762136, 0.762137], [1.5864, 1.58641], [-1.50617, -1.50616]]
+  560.079 μs (4832 allocations: 199.20 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-2.79797, -2.79796], [-1.73355, -1.73354], [-0.606508, -0.606507], [-2.44883, -2.44882], [1.71346, 1.71347], [-2.72265, -2.72264], [-1.55253, -1.55252], [0.289974, 0.289975], [0.252046, 0.252047], [0.762136, 0.762137], [1.5864, 1.58641], [-1.50617, -1.50616]]
+
+For n = 13
+  656.259 μs (5226 allocations: 218.75 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[2.82409, 2.8241], [-1.48454, -1.48453], [-0.31499, -0.314989], [3.10565, 3.10566], [1.40989, 1.4099], [1.96286, 1.96287], [1.07219, 1.0722], [-4.11245, -4.11244], [-1.54954, -1.54953], [1.06494, 1.06495], [-3.6406, -3.64059], [1.50881, 1.50882], [1.19286, 1.19287]]
+  676.472 μs (5232 allocations: 216.09 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[2.82409, 2.8241], [-1.48454, -1.48453], [-0.31499, -0.314989], [3.10565, 3.10566], [1.40989, 1.4099], [1.96286, 1.96287], [1.07219, 1.0722], [-4.11245, -4.11244], [-1.54954, -1.54953], [1.06494, 1.06495], [-3.6406, -3.64059], [1.50881, 1.50882], [1.19286, 1.19287]]
+
+For n = 14
+  780.545 μs (5626 allocations: 236.19 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[0.78063, 0.780631], [0.0801395, 0.0801396], [-0.619735, -0.619734], [-0.408535, -0.408534], [-0.35918, -0.359179], [-0.279269, -0.279268], [-0.381064, -0.381063], [-0.583649, -0.583648], [0.15355, 0.153551], [0.303571, 0.303572], [-0.263495, -0.263494], [0.745922, 0.745923], [-1.25136, -1.25135], [-0.306904, -0.306903]]
+  792.250 μs (5632 allocations: 233.17 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[0.78063, 0.780631], [0.0801395, 0.0801396], [-0.619735, -0.619734], [-0.408535, -0.408534], [-0.35918, -0.359179], [-0.279269, -0.279268], [-0.381064, -0.381063], [-0.583649, -0.583648], [0.15355, 0.153551], [0.303571, 0.303572], [-0.263495, -0.263494], [0.745922, 0.745923], [-1.25136, -1.25135], [-0.306904, -0.306903]]
+
+For n = 15
+  896.447 μs (6026 allocations: 253.50 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-1.31651, -1.3165], [1.50793, 1.50794], [-1.66271, -1.6627], [0.42862, 0.428621], [0.764462, 0.764463], [0.376219, 0.37622], [0.0408481, 0.0408482], [0.078441, 0.0784411], [-1.53466, -1.53465], [0.498318, 0.498319], [-0.257995, -0.257994], [-0.604852, -0.604851], [-0.386613, -0.386612], [-0.0438574, -0.0438573], [0.440631, 0.440632]]
+  950.789 μs (6032 allocations: 250.02 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-1.31651, -1.3165], [1.50793, 1.50794], [-1.66271, -1.6627], [0.42862, 0.428621], [0.764462, 0.764463], [0.376219, 0.37622], [0.0408481, 0.0408482], [0.078441, 0.0784411], [-1.53466, -1.53465], [0.498318, 0.498319], [-0.257995, -0.257994], [-0.604852, -0.604851], [-0.386613, -0.386612], [-0.0438574, -0.0438573], [0.440631, 0.440632]]
+
+For n = 16
+  1.076 ms (6426 allocations: 270.61 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-0.0621213, -0.0621212], [-0.436339, -0.436338], [0.680463, 0.680464], [0.55125, 0.551251], [0.979387, 0.979388], [-0.904311, -0.90431], [-1.34208, -1.34207], [2.02843, 2.02844], [1.99408, 1.99409], [-1.1726, -1.17259], [1.25116, 1.25117], [-0.313683, -0.313682], [1.09734, 1.09735], [0.232912, 0.232913], [-1.68377, -1.68376], [-0.425297, -0.425296]]
+  1.116 ms (6432 allocations: 266.64 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-0.0621213, -0.0621212], [-0.436339, -0.436338], [0.680463, 0.680464], [0.55125, 0.551251], [0.979387, 0.979388], [-0.904311, -0.90431], [-1.34208, -1.34207], [2.02843, 2.02844], [1.99408, 1.99409], [-1.1726, -1.17259], [1.25116, 1.25117], [-0.313683, -0.313682], [1.09734, 1.09735], [0.232912, 0.232913], [-1.68377, -1.68376], [-0.425297, -0.425296]]
+
+For n = 17
+  1.205 ms (6826 allocations: 288.27 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[0.543059, 0.54306], [-1.29556, -1.29555], [-1.7007, -1.70069], [1.43802, 1.43803], [0.179184, 0.179185], [1.69012, 1.69013], [1.34535, 1.34536], [1.49156, 1.49157], [1.15319, 1.1532], [-0.629132, -0.629131], [-1.15364, -1.15363], [0.627152, 0.627153], [-1.82, -1.81999], [-2.69038, -2.69037], [2.13837, 2.13838], [-0.653879, -0.653878], [-0.585587, -0.585586]]
+  1.284 ms (6832 allocations: 283.75 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[0.543059, 0.54306], [-1.29556, -1.29555], [-1.7007, -1.70069], [1.43802, 1.43803], [0.179184, 0.179185], [1.69012, 1.69013], [1.34535, 1.34536], [1.49156, 1.49157], [1.15319, 1.1532], [-0.629132, -0.629131], [-1.15364, -1.15363], [0.627152, 0.627153], [-1.82, -1.81999], [-2.69038, -2.69037], [2.13837, 2.13838], [-0.653879, -0.653878], [-0.585587, -0.585586]]
+
+For n = 18
+  1.541 ms (7226 allocations: 305.64 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-4.4038, -4.40379], [35.4039, 35.404], [23.2018, 23.2019], [29.5486, 29.5487], [-4.1716, -4.17159], [-10.229, -10.2289], [-4.36794, -4.36793], [-3.06734, -3.06733], [28.4815, 28.4816], [20.6252, 20.6253], [8.52634, 8.52635], [19.3466, 19.3467], [-17.4713, -17.4712], [-40.2325, -40.2324], [8.11921, 8.11922], [-23.8136, -23.8135], [17.8947, 17.8948], [55.5122, 55.5123]]
+  1.666 ms (7232 allocations: 300.63 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-4.4038, -4.40379], [35.4039, 35.404], [23.2018, 23.2019], [29.5486, 29.5487], [-4.1716, -4.17159], [-10.229, -10.2289], [-4.36794, -4.36793], [-3.06734, -3.06733], [28.4815, 28.4816], [20.6252, 20.6253], [8.52634, 8.52635], [19.3466, 19.3467], [-17.4713, -17.4712], [-40.2325, -40.2324], [8.11921, 8.11922], [-23.8136, -23.8135], [17.8947, 17.8948], [55.5122, 55.5123]]
+
+For n = 19
+  1.769 ms (7626 allocations: 323.36 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-0.16665, -0.166649], [-0.836062, -0.836061], [0.191566, 0.191567], [0.359444, 0.359445], [0.489007, 0.489008], [-0.125, -0.124999], [-0.397535, -0.397534], [-0.45875, -0.458749], [0.735181, 0.735182], [1.68832, 1.68833], [0.555231, 0.555232], [1.00734, 1.00735], [-0.433021, -0.43302], [0.742541, 0.742542], [-1.01507, -1.01506], [1.42938, 1.42939], [-0.306669, -0.306668], [-0.657369, -0.657368], [-0.139753, -0.139752]]
+  1.930 ms (7632 allocations: 317.73 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-0.16665, -0.166649], [-0.836062, -0.836061], [0.191566, 0.191567], [0.359444, 0.359445], [0.489007, 0.489008], [-0.125, -0.124999], [-0.397535, -0.397534], [-0.45875, -0.458749], [0.735181, 0.735182], [1.68832, 1.68833], [0.555231, 0.555232], [1.00734, 1.00735], [-0.433021, -0.43302], [0.742541, 0.742542], [-1.01507, -1.01506], [1.42938, 1.42939], [-0.306669, -0.306668], [-0.657369, -0.657368], [-0.139753, -0.139752]]
+
+For n = 20
+  1.964 ms (8026 allocations: 340.89 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[1.35974, 1.35975], [0.0315064, 0.0315065], [1.62038, 1.62039], [1.49915, 1.49916], [0.846944, 0.846945], [1.45553, 1.45554], [4.53343, 4.53344], [1.90342, 1.90343], [-4.99214, -4.99213], [2.07785, 2.07786], [-0.444216, -0.444215], [0.736985, 0.736986], [1.07722, 1.07723], [2.75322, 2.75323], [1.59534, 1.59535], [2.78334, 2.78335], [-5.37279, -5.37278], [-0.843305, -0.843304], [2.58871, 2.58872], [0.805166, 0.805167]]
+  2.225 ms (8032 allocations: 334.67 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[1.35974, 1.35975], [0.0315064, 0.0315065], [1.62038, 1.62039], [1.49915, 1.49916], [0.846944, 0.846945], [1.45553, 1.45554], [4.53343, 4.53344], [1.90342, 1.90343], [-4.99214, -4.99213], [2.07785, 2.07786], [-0.444216, -0.444215], [0.736985, 0.736986], [1.07722, 1.07723], [2.75322, 2.75323], [1.59534, 1.59535], [2.78334, 2.78335], [-5.37279, -5.37278], [-0.843305, -0.843304], [2.58871, 2.58872], [0.805166, 0.805167]]
+
+For n = 21
+  2.279 ms (8426 allocations: 358.86 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[0.531469, 0.53147], [-1.55322, -1.55321], [0.610774, 0.610775], [-0.972668, -0.972667], [-0.995764, -0.995763], [1.53137, 1.53138], [0.226348, 0.226349], [-0.0360293, -0.0360292], [0.300785, 0.300786], [1.16974, 1.16975], [0.224403, 0.224404], [-1.59163, -1.59162], [-0.33666, -0.336659], [-0.365669, -0.365668], [0.671668, 0.671669], [0.166757, 0.166758], [-0.0187087, -0.0187086], [0.000503947, 0.000503948], [0.747613, 0.747614], [-0.436936, -0.436935], [-1.07276, -1.07275]]
+  2.523 ms (8432 allocations: 351.97 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[0.531469, 0.53147], [-1.55322, -1.55321], [0.610774, 0.610775], [-0.972668, -0.972667], [-0.995764, -0.995763], [1.53137, 1.53138], [0.226348, 0.226349], [-0.0360293, -0.0360292], [0.300785, 0.300786], [1.16974, 1.16975], [0.224403, 0.224404], [-1.59163, -1.59162], [-0.33666, -0.336659], [-0.365669, -0.365668], [0.671668, 0.671669], [0.166757, 0.166758], [-0.0187087, -0.0187086], [0.000503947, 0.000503948], [0.747613, 0.747614], [-0.436936, -0.436935], [-1.07276, -1.07275]]
+
+For n = 22
+  2.635 ms (8826 allocations: 376.55 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-0.0576031, -0.057603], [-0.19045, -0.190449], [-0.145126, -0.145125], [0.0613842, 0.0613843], [-0.0716046, -0.0716045], [-0.020995, -0.0209949], [-0.399329, -0.399328], [-0.291664, -0.291663], [-0.131188, -0.131187], [0.287153, 0.287154], [-0.476477, -0.476476], [0.937626, 0.937627], [-0.38509, -0.385089], [0.0836971, 0.0836972], [0.176637, 0.176638], [0.887916, 0.887917], [0.0787311, 0.0787312], [-0.404081, -0.40408], [-0.319168, -0.319167], [-0.981298, -0.981297], [-0.0197936, -0.0197935], [0.589842, 0.589843]]
+  2.997 ms (8832 allocations: 369.03 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-0.0576031, -0.057603], [-0.19045, -0.190449], [-0.145126, -0.145125], [0.0613842, 0.0613843], [-0.0716046, -0.0716045], [-0.020995, -0.0209949], [-0.399329, -0.399328], [-0.291664, -0.291663], [-0.131188, -0.131187], [0.287153, 0.287154], [-0.476477, -0.476476], [0.937626, 0.937627], [-0.38509, -0.385089], [0.0836971, 0.0836972], [0.176637, 0.176638], [0.887916, 0.887917], [0.0787311, 0.0787312], [-0.404081, -0.40408], [-0.319168, -0.319167], [-0.981298, -0.981297], [-0.0197936, -0.0197935], [0.589842, 0.589843]]
+
+For n = 23
+  2.810 ms (9226 allocations: 394.70 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-1.89769, -1.89768], [-2.13486, -2.13485], [-0.0239536, -0.0239535], [2.26537, 2.26538], [0.142881, 0.142882], [-1.55981, -1.5598], [-1.5167, -1.51669], [-0.511803, -0.511802], [0.571971, 0.571972], [-0.00964798, -0.00964797], [-0.858116, -0.858115], [-0.239622, -0.239621], [-0.891741, -0.89174], [-2.30866, -2.30865], [-1.76962, -1.76961], [-2.61382, -2.61381], [-0.457754, -0.457753], [-1.03219, -1.03218], [-1.18658, -1.18657], [0.0264487, 0.0264488], [0.0913128, 0.0913129],[-0.0189498, -0.0189497], [-0.922417, -0.922416]]
+  3.222 ms (9232 allocations: 386.45 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-1.89769, -1.89768], [-2.13486, -2.13485], [-0.0239536, -0.0239535], [2.26537, 2.26538], [0.142881, 0.142882], [-1.55981, -1.5598], [-1.5167, -1.51669], [-0.511803, -0.511802], [0.571971, 0.571972], [-0.00964798, -0.00964797], [-0.858116, -0.858115], [-0.239622, -0.239621], [-0.891741, -0.89174], [-2.30866, -2.30865], [-1.76962, -1.76961], [-2.61382, -2.61381], [-0.457754, -0.457753], [-1.03219, -1.03218], [-1.18658, -1.18657], [0.0264487, 0.0264488], [0.0913128, 0.0913129], [-0.0189498, -0.0189497], [-0.922417, -0.922416]]
+
+For n = 24
+  3.219 ms (9626 allocations: 412.64 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-0.124356, -0.124355], [-0.00958074, -0.00958073], [0.539476, 0.539477], [-0.388642, -0.388641], [0.625504, 0.625505], [1.09303, 1.09304], [-1.11944, -1.11943], [0.638843, 0.638844], [-0.522925, -0.522924], [-1.37454, -1.37453], [0.642279, 0.64228], [0.764112, 0.764113], [-0.568954, -0.568953], [-0.558264, -0.558263], [0.58252, 0.582521], [0.735529, 0.73553], [-0.365248, -0.365247], [-1.32431, -1.3243], [-0.150569, -0.150568], [-0.399501, -0.3995], [1.37909, 1.3791], [0.00529648, 0.00529649], [0.21186, 0.211861], [0.6825, 0.682501]]
+  3.707 ms (9632 allocations: 403.56 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-0.124356, -0.124355], [-0.00958074, -0.00958073], [0.539476, 0.539477], [-0.388642, -0.388641], [0.625504, 0.625505], [1.09303,1.09304], [-1.11944, -1.11943], [0.638843, 0.638844], [-0.522925, -0.522924], [-1.37454, -1.37453], [0.642279, 0.64228], [0.764112, 0.764113], [-0.568954, -0.568953], [-0.558264, -0.558263], [0.58252, 0.582521], [0.735529, 0.73553], [-0.365248, -0.365247], [-1.32431, -1.3243], [-0.150569, -0.150568], [-0.399501, -0.3995], [1.37909, 1.3791], [0.00529648, 0.00529649], [0.21186, 0.211861], [0.6825, 0.682501]]
+
+For n = 25
+  3.587 ms (10028 allocations: 431.00 KiB)
+Array: IntervalArithmetic.Interval{Float64}[[-1.74061, -1.7406], [0.412884, 0.412885], [0.623495, 0.623496], [-1.42243, -1.42242], [-0.774224, -0.774223], [-0.0379006, -0.0379005], [0.735943, 0.735944], [0.0682942, 0.0682943], [0.992878, 0.992879], [-0.492704, -0.492703], [-1.38519, -1.38518], [0.387523, 0.387524], [-3.107, -3.10699], [0.0425826, 0.0425827], [-1.00909, -1.00908], [-0.164128, -0.164127], [0.0325966, 0.0325967], [-0.52808, -0.528079], [0.499164, 0.499165], [-0.00699219, -0.00699218], [0.034, 0.0340001], [0.878137, 0.878138], [0.396617, 0.396618], [-1.18615, -1.18614], [0.043228, 0.0432281]]
+  4.199 ms (10032 allocations: 421.02 KiB)
+MArray: IntervalArithmetic.Interval{Float64}[[-1.74061, -1.7406], [0.412884, 0.412885], [0.623495, 0.623496], [-1.42243, -1.42242], [-0.774224, -0.774223], [-0.0379006, -0.0379005], [0.735943, 0.735944], [0.0682942, 0.0682943], [0.992878, 0.992879], [-0.492704, -0.492703], [-1.38519, -1.38518], [0.387523, 0.387524], [-3.107, -3.10699], [0.0425826,0.0425827], [-1.00909, -1.00908], [-0.164128, -0.164127], [0.0325966, 0.0325967], [-0.52808, -0.528079], [0.499164, 0.499165], [-0.00699219, -0.00699218], [0.034, 0.0340001],[0.878137, 0.878138], [0.396617, 0.396618], [-1.18615, -1.18614], [0.043228, 0.0432281]]
diff --git a/perf/linear_eq_tabular.txt b/perf/linear_eq_tabular.txt
new file mode 100644
index 00000000..69028385
--- /dev/null
+++ b/perf/linear_eq_tabular.txt
@@ -0,0 +1,49 @@
+Gauss Seidel
+
+│ Row │ variable │ Array       │ Contractor  │ MArray      │ SArray1     │ SArray2     │
+├─────┼──────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┤
+│ 1   │ n = 1    │ 1.1426e-5   │ 2.4173e-5   │ 1.5803e-5   │ 1.014e-5    │ 9.512e-6    │
+│ 2   │ n = 2    │ 2.4413e-5   │ 4.0504e-5   │ 3.9829e-5   │ 2.4083e-5   │ 2.4083e-5   │
+│ 3   │ n = 3    │ 4.5194e-5   │ 6.4343e-5   │ 9.1236e-5   │ 4.5533e-5   │ 5.1221e-5   │
+│ 4   │ n = 4    │ 7.2639e-5   │ 9.3671e-5   │ 0.00016058  │ 7.0972e-5   │ 6.9937e-5   │
+│ 5   │ n = 5    │ 0.000103281 │ 0.000127489 │ 0.000262814 │ 0.000104823 │ 0.000102729 │
+│ 6   │ n = 6    │ 0.000141315 │ 0.000169992 │ 0.000416266 │ 0.000144021 │ 0.000139086 │
+│ 7   │ n = 7    │ 0.000195001 │ 0.000217059 │ 0.000680848 │ 0.000198164 │ 0.000191145 │
+│ 8   │ n = 8    │ 0.000247193 │ 0.00027535  │ 0.00122559  │ 0.000251089 │ 0.000241972 │
+│ 9   │ n = 9    │ 0.000306262 │ 0.000338028 │ 0.0020223   │ 0.000310961 │ 0.000299153 │
+│ 10  │ n = 10   │ 0.00037781  │ 0.000414073 │ 0.00294134  │ 0.000386235 │ 0.000367335 │
+│ 11  │ n = 11   │ 0.000465017 │ 0.000489114 │ 0.00585246  │ 0.000473425 │ 0.000447176 │
+│ 12  │ n = 12   │ 0.00055403  │ 0.000587152 │ 0.0110894   │ 0.000573395 │ 0.00054274  │
+│ 13  │ n = 13   │ 0.000659994 │ 0.000692402 │ 0.0160662   │ 0.000686098 │ 0.000643391 │
+│ 14  │ n = 14   │ 0.000777718 │ 0.000793626 │ 0.0184595   │ 0.000800443 │ 0.000767576 │
+│ 15  │ n = 15   │ 0.000910174 │ 0.000952207 │ 0.0239422   │ 0.000972855 │ 0.000899402 │
+│ 16  │ n = 16   │ 0.00107785  │ 0.00111443  │ 0.0303295   │ 0.00113467  │ 0.0010867   │
+│ 17  │ n = 17   │ 0.00122326  │ 0.00125759  │ 0.0398127   │ 0.00134961  │ 0.00125675  │
+│ 18  │ n = 18   │ 0.00167176  │ 0.00159463  │ 0.043795    │ 0.00181409  │ 0.00162584  │
+│ 19  │ n = 19   │ 0.0018632   │ 0.00185617  │ 0.0549107   │ 0.00208529  │ 0.00196587  │
+│ 20  │ n = 20   │ 0.0020604   │ 0.00210595  │ 0.0635598   │ 0.00232704  │ 0.00212374  │
+
+10×7 DataFrames.DataFrame
+│ Row │ n  │ Array       │ Contractor  │ ContractorMArray │ ContractorSArray │ MArray      │ SArray      │
+├─────┼────┼─────────────┼─────────────┼──────────────────┼──────────────────┼─────────────┼─────────────┤
+│ 1   │ 1  │ 1.0188e-5   │ 2.5276e-5   │ 0.000781355      │ 0.000764627      │ 1.6959e-5   │ 1.1216e-5   │
+│ 2   │ 10 │ 0.000377853 │ 0.000407817 │ 0.00122723       │ 0.00123259       │ 0.0029247   │ 0.000393896 │
+│ 3   │ 2  │ 2.609e-5    │ 4.3597e-5   │ 0.000802467      │ 0.000796962      │ 4.1039e-5   │ 2.6654e-5   │
+│ 4   │ 3  │ 4.7362e-5   │ 6.317e-5    │ 0.000832029      │ 0.000821289      │ 9.1695e-5   │ 4.7007e-5   │
+│ 5   │ 4  │ 7.0378e-5   │ 9.1502e-5   │ 0.00085429       │ 0.000844048      │ 0.000163604 │ 7.1166e-5   │
+│ 6   │ 5  │ 0.000104813 │ 0.000129112 │ 0.000898626      │ 0.000889969      │ 0.000268756 │ 0.000112476 │
+│ 7   │ 6  │ 0.000142316 │ 0.000168891 │ 0.000941422      │ 0.000922727      │ 0.000430442 │ 0.000146906 │
+│ 8   │ 7  │ 0.000192073 │ 0.000218861 │ 0.00108794       │ 0.000997821      │ 0.00073304  │ 0.000198962 │
+│ 9   │ 8  │ 0.000242832 │ 0.000273341 │ 0.00107159       │ 0.00106769       │ 0.00126407  │ 0.000250956 │
+│ 10  │ 9  │ 0.000308256 │ 0.000333722 │ 0.00116203       │ 0.00116428       │ 0.00196451  │ 0.000314499 │
+
+Gauss Elimination
+
+5×3 DataFrames.DataFrame
+│ Row │ n │ Base.\\    │ Gauss Elimination │
+├─────┼───┼────────────┼───────────────────┤
+│ 1   │ 1 │ 1.84e-6    │ 8.127e-6          │
+│ 2   │ 2 │ 3.1615e-6  │ 1.0029e-5         │
+│ 3   │ 3 │ 4.53986e-6 │ 1.1211e-5         │
+│ 4   │ 4 │ 6.8655e-6  │ 1.4342e-5         │
+│ 5   │ 5 │ 1.0773e-5  │ 1.7725e-5         │
diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl
index a363b35c..480a4345 100644
--- a/src/IntervalRootFinding.jl
+++ b/src/IntervalRootFinding.jl
@@ -9,7 +9,7 @@ using ForwardDiff
 using StaticArrays
 
 
-import Base: ⊆, show, big
+import Base: ⊆, show, big, \
 
 import Polynomials: roots
 
@@ -18,7 +18,10 @@ export
     derivative, jacobian,  # reexport derivative from ForwardDiff
     Root, is_unique,
     roots, find_roots,
-    bisect, newton1d, quadratic_roots
+    bisect, newton1d, quadratic_roots,
+    gauss_seidel_interval, gauss_seidel_interval!,
+    gauss_seidel_contractor, gauss_seidel_contractor!,
+    gauss_elimination_interval, gauss_elimination_interval!
 
 export isunique, root_status
 
@@ -59,6 +62,7 @@ include("contractors.jl")
 include("roots.jl")
 include("newton1d.jl")
 include("quadratic.jl")
+include("linear_eq.jl")
 
 
 gradient(f) = X -> ForwardDiff.gradient(f, SVector(X))
diff --git a/src/linear_eq.jl b/src/linear_eq.jl
new file mode 100644
index 00000000..b4242b17
--- /dev/null
+++ b/src/linear_eq.jl
@@ -0,0 +1,160 @@
+"""
+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)
+
+    n = size(A, 1)
+    x = similar(b)
+    x .= -1e16..1e16
+    gauss_seidel_interval!(x, A, b, precondition=precondition, maxiter=maxiter)
+    return x
+end
+"""
+Iteratively solves the system of interval linear
+equations and returns the solution set. Uses the
+Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
+Keyword `precondition` to turn preconditioning off.
+Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
+"""
+function gauss_seidel_interval!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+
+    precondition && ((A, b) = preconditioner(A, b))
+
+    n = size(A, 1)
+
+    @inbounds for iter in 1:maxiter
+        x¹ = copy(x)
+        for i in 1:n
+            Y = b[i]
+            for j in 1:n
+                (i == j) || (Y -= A[i, j] * x[j])
+            end
+            Z = extended_div(Y, A[i, i])
+            x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2])
+        end
+        if all(x .== x¹)
+            break
+        end
+    end
+    x
+end
+
+function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+
+    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)
+
+    precondition && ((A, b) = preconditioner(A, b))
+
+    n = size(A, 1)
+
+    diagA = Diagonal(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
+    end
+    inv_diagA = inv(diagA)
+
+    for iter in 1:maxiter
+        x¹ = copy(x)
+        x = x .∩ (inv_diagA * (b - extdiagA * x))
+        if all(x .== x¹)
+            break
+        end
+    end
+    x
+end
+
+function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true)
+
+    n = size(A, 1)
+    x = similar(b)
+    x .= -1e16..1e16
+    x = gauss_elimination_interval!(x, A, b, precondition=precondition)
+
+    return x
+end
+"""
+Solves the system of linear equations using Gaussian Elimination,
+with (or without) preconditioning. (kwarg - `precondition`)
+Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis - Page 72
+"""
+function gauss_elimination_interval!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true)
+
+    if precondition
+        ((a, b) = preconditioner(a, b))
+    else
+        a = copy(a)
+        b = copy(b)
+    end
+    n = size(a, 1)
+
+    p = zeros(x)
+
+    for i in 1:(n-1)
+        if 0 ∈ a[i, i] # diagonal matrix is not invertible
+            p .= entireinterval(b[1])
+            return p .∩ x
+        end
+
+        for j in (i+1):n
+            α = a[j, i] / a[i, i]
+            b[j] -= α * b[i]
+
+            for k in (i+1):n
+                a[j, k] -= α * a[i, k]
+            end
+        end
+    end
+
+    for i in n:-1:1
+        sum = 0
+        for j in (i+1):n
+            sum += a[i, j] * p[j]
+        end
+        p[i] = (b[i] - sum) / a[i, i]
+    end
+
+    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...)
diff --git a/test/linear_eq.jl b/test/linear_eq.jl
new file mode 100644
index 00000000..2cc49d08
--- /dev/null
+++ b/test/linear_eq.jl
@@ -0,0 +1,40 @@
+using IntervalArithmetic, StaticArrays, IntervalRootFinding
+
+function randVec(n::Int)
+    a = randn(n)
+    A = Interval.(a)
+    mA = MVector{n}(A)
+    sA = SVector{n}(A)
+    return A, mA, sA
+end
+
+function randMat(n::Int)
+    a = randn(n, n)
+    A = Interval.(a)
+    mA = MMatrix{n, n}(A)
+    sA = SMatrix{n, n}(A)
+    return A, mA, sA
+end
+
+@testset "Linear Equations" begin
+
+    A = [[2..3 0..1;1..2 2..3], ]
+    b = [[0..120, 60..240], ]
+    x = [[-120..90, -60..240], ]
+
+    for i in 1:10
+        rand_a = randMat(i)[1]
+        rand_b = randVec(i)[1]
+        rand_c = rand_a * rand_b
+        push!(A, rand_a)
+        push!(b, rand_c)
+        push!(x, rand_b)
+    end
+    for i in 1:length(A)
+        for precondition in (false, true)
+            for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval, \)
+                @test all(x[i] .⊆ f(A[i], b[i]))
+            end
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index bf7d49a3..e514cec7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,3 +7,4 @@ include("roots.jl")
 include("newton1d.jl")
 include("quadratic.jl")
 include("test_smiley.jl")
+include("linear_eq.jl")