diff --git a/contents/gaussian_elimination/code/haskell/gaussianElimination.hs b/contents/gaussian_elimination/code/haskell/gaussianElimination.hs
new file mode 100644
index 000000000..62da17926
--- /dev/null
+++ b/contents/gaussian_elimination/code/haskell/gaussianElimination.hs
@@ -0,0 +1,89 @@
+import Data.Array
+import Data.Function (on)
+import Data.List (intercalate, maximumBy)
+import Data.Ratio
+
+type Matrix a = Array (Int, Int) a
+
+type Vector a = Array Int a
+
+swapRows :: Int -> Int -> Matrix a -> Matrix a
+swapRows r1 r2 m
+  | r1 == r2 = m
+  | otherwise =
+    m //
+    concat [[((r2, c), m ! (r1, c)), ((r1, c), m ! (r2, c))] | c <- [c1 .. cn]]
+  where
+    ((_, c1), (_, cn)) = bounds m
+
+subRows ::
+     Fractional a
+  => (Int, Int) -- pivot location
+  -> (Int, Int) -- rows to cover
+  -> (Int, Int) -- columns to cover
+  -> Matrix a
+  -> Matrix a
+subRows (r, c) (r1, rn) (c1, cn) m =
+  accum
+    (-)
+    m
+    [ ((i, j), m ! (i, c) * m ! (r, j) / m ! (r, c))
+    | i <- [r1 .. rn]
+    , j <- [c1 .. cn]
+    ]
+
+gaussianElimination :: (Fractional a, Ord a) => Matrix a -> Matrix a
+gaussianElimination mat = go (r1, c1) mat
+  where
+    ((r1, c1), (rn, cn)) = bounds mat
+    go (r, c) m
+      | c == cn = m
+      | pivot == 0 = go (r, c + 1) m
+      | otherwise = go (r + 1, c + 1) $ subRows (r, c) (r + 1, rn) (c, cn) m'
+      where
+        (target, pivot) =
+          maximumBy (compare `on` abs . snd) [(k, m ! (k, c)) | k <- [r .. rn]]
+        m' = swapRows r target m
+
+gaussJordan :: (Fractional a, Eq a) => Matrix a -> Matrix a
+gaussJordan mat = go (r1, c1) mat
+  where
+    ((r1, c1), (rn, cn)) = bounds mat
+    go (r, c) m
+      | c == cn = m
+      | m ! (r, c) == 0 = go (r, c + 1) m
+      | otherwise = go (r + 1, c + 1) $ subRows (r, c) (r1, r - 1) (c, cn) m'
+      where
+        m' = accum (/) m [((r, j), m ! (r, c)) | j <- [c .. cn]]
+
+backSubstitution :: (Fractional a) => Matrix a -> Vector a
+backSubstitution m = sol
+  where
+    ((r1, _), (rn, cn)) = bounds m
+    sol =
+      listArray (r1, rn) [(m ! (r, cn) - sum' r) / m ! (r, r) | r <- [r1 .. rn]]
+    sum' r = sum [m ! (r, k) * sol ! k | k <- [r + 1 .. rn]]
+
+printM :: (Show a) => Matrix a -> String
+printM m =
+  let ((r1, c1), (rn, cn)) = bounds m
+   in unlines
+        [ intercalate "\t" [show $ m ! (r, c) | c <- [c1 .. cn]]
+        | r <- [r1 .. rn]
+        ]
+
+printV :: (Show a) => Vector a -> String
+printV = unlines . map show . elems
+
+main :: IO ()
+main = do
+  let mat = [2, 3, 4, 6, 1, 2, 3, 4, 3, -4, 0, 10] :: [Ratio Int]
+      m = listArray ((1, 1), (3, 4)) mat
+  putStrLn "Original Matrix:"
+  putStrLn $ printM m
+  putStrLn "Echelon form"
+  putStrLn $ printM $ gaussianElimination m
+  putStrLn "Reduced echelon form"
+  putStrLn $ printM $ gaussJordan $ gaussianElimination m
+  putStrLn "Solution from back substitution"
+  putStrLn $ printV $ backSubstitution $ gaussianElimination m
diff --git a/contents/gaussian_elimination/gaussian_elimination.md b/contents/gaussian_elimination/gaussian_elimination.md
index 6548edd05..f7e31bfdb 100644
--- a/contents/gaussian_elimination/gaussian_elimination.md
+++ b/contents/gaussian_elimination/gaussian_elimination.md
@@ -149,7 +149,7 @@ $$
 \begin{array}{ccc|c}
 1 & 0 & 0 & \frac{18}{11} \\
 0 & 1 & 0 & \frac{-14}{11} \\
-0 & 0 & 1 & \frac{18}{11} 
+0 & 0 & 1 & \frac{18}{11}
 \end{array}
 \right]
 $$
@@ -358,6 +358,8 @@ In code, this looks like:
 [import:13-44, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs)
+{% sample lang="hs" %}
+[import:10-36, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 Now, to be clear: this algorithm creates an upper-triangular matrix.
@@ -390,6 +392,8 @@ This code does not exist yet in C, so here's Julia code (sorry for the inconveni
 {% sample lang="rs" %}
 This code does not exist yet in rust, so here's Julia code (sorry for the inconvenience)
 [import:70-96, lang:"julia"](code/julia/gaussian_elimination.jl)
+{% sample lang="hs" %}
+[import:38-46, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 ## Back-substitution
@@ -418,6 +422,8 @@ In code, this involves keeping a rolling sum of all the values we substitute in
 [import:46-58, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs)
+{% sample lang="hs" %}
+[import:48-53, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 ## Conclusions
@@ -438,6 +444,8 @@ As for what's next... Well, we are in for a treat! The above algorithm clearly h
 [import, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import, lang:"rust"](code/rust/gaussian_elimination.rs)
+{% sample lang="hs" %}
+[import, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 
@@ -463,4 +471,3 @@ $$
 \newcommand{\bfomega}{\boldsymbol{\omega}}
 \newcommand{\bftau}{\boldsymbol{\tau}}
 $$
-