diff --git a/sm/mat.cppm b/sm/mat.cppm
index e8cedcf..b0e20b8 100644
--- a/sm/mat.cppm
+++ b/sm/mat.cppm
@@ -869,8 +869,7 @@ export namespace sm
 
             for (std::uint32_t k = 1u; k <= Nr; ++k) {
                 M = (*this) * M;  // M_k = A * M_{k-1}, where M_0 = I
-                F trace = M.trace();
-                coeffs[Nr - k] = -trace / F(k);
+                coeffs[Nr - k] = -M.trace() / static_cast<F>(k);
 
                 if (k < Nr) {
                     // M = M_k + c_{n-k} * I
@@ -984,52 +983,45 @@ export namespace sm
             return rref;
         }
 
-        // Solve matrix by back-substitution. row_echelon_form is fine (don't need reduced row
-        // echelon form). Nans in result mean there was no solution.
+        /*!
+         * Using back-substitution, solve the equation Ax = B, where A and B are represented in
+         * *this in the usual way, with A being a square matrix on the left and B the right-most
+         * column. *this should thus be the 'augmented matrix'. This matrix should have been reduced
+         * to row_echelon_form before calling back_substitution().
+         *
+         * \return A vec of the solutions, x. NaNs in result mean there was no solution. Currently,
+         * I don't represent a free variable. Is that possible? I *could* use limits::max?
+         */
         template<typename Fy=F> requires (Nc == Nr + 1)
-        sm::vec<F, Nr> back_substitution()
+        sm::vec<F, Nr> back_substitution() const noexcept
         {
             using F_el = typename value_type_of<F>::type;
             constexpr F_el eps = F_el{1e-14}; // needs to be F_el if F is complex
 
             sm::vec<F, Nr> x = {};
-            constexpr bool human_readable_version = false;
-
-            if constexpr (human_readable_version) {
-                // Developed first with human-friendly matrix array indexing like A(r, c)
-                for (std::uint32_t i = (Nr - 1u); i != std::numeric_limits<std::uint32_t>::max(); i--) {
-                    if (std::abs ((*this)(i, i)) > eps) {   // abs value of the triangular element
-                        const std::uint32_t n = Nr - i - 1; // number of elements in the sum
-                        F sum = F{0};
-                        for (std::uint32_t j = 0; j < n; ++j) {
-                            sum += (*this)(i, i + n - j) * x[i + n - j];
-                        }
-                        x[i] = ((*this)(i, Nc - 1) - sum ) / (*this)(i, i);
-                    }  // else x[i] remains 0. This could mean no solution if (*this)(i, Nc - 1) != 0
-                }
-            } else {
-                // But compile this; Converting operator indexing to array indices as: (r, c) -> idx
-                // = (c * Nr) + r and removing need for the sum variable.
-                for (std::uint32_t i = (Nr - 1u); i != std::numeric_limits<std::uint32_t>::max(); i--) {
-                    const std::uint32_t de = (i * Nr) + i;        // diagonal element
-                    const std::uint32_t le = ((Nc - 1) * Nr) + i; // last element
-                    if (std::abs (this->arr[de]) > eps) {
-                        const std::uint32_t n = Nr - i - 1;
-                        x[i] = this->arr[le];
-                        for (std::uint32_t j = 0; j < n; ++j) {
-                            const std::uint32_t se = (Nr - 1 - j) * Nr + i;
-                            x[i] -= this->arr[se] * x[i + n - j];
-                        }
-                        x[i] /= this->arr[de];
+
+            for (std::uint32_t i = (Nr - 1u); i != std::numeric_limits<std::uint32_t>::max(); i--) {
+                const std::uint32_t de = (i * Nr) + i;        // diagonal element A(r, r)
+                const std::uint32_t le = ((Nc - 1) * Nr) + i; // last element A(r, Nc - 1)
+                if (std::abs (this->arr[de]) > eps) { // std::abs not sm::cem::abs as F may be complex
+                    const std::uint32_t n = Nr - i - 1;
+                    x[i] = this->arr[le];
+                    for (std::uint32_t j = 0; j < n; ++j) {
+                        x[i] -= this->arr[(Nr - 1 - j) * Nr + i] * x[i + n - j]; // sum element is A(i, i + n - j)
+                    }
+                    x[i] /= this->arr[de];
+                } else {
+                    // arr[de] == 0. If the value in the right-most column is 0 AND the value in the
+                    // diagonal is 0, then we have a free value and x[i] remains 0. If the value in
+                    // the right most column is != 0, but arr[de] == 0, we have no solution.
+                    if constexpr (sm::is_complex<F>::value) {
+                        if (this->arr[le] != F{0}) { x[i] = F{std::numeric_limits<F_el>::quiet_NaN()}; }
                     } else {
-                        if constexpr (sm::is_complex<F>::value) {
-                            if (this->arr[le] != F{0}) { x[i] = F{std::numeric_limits<F_el>::quiet_NaN()}; }
-                        } else {
-                            if (this->arr[le] != F{0}) { x[i] = std::numeric_limits<F>::quiet_NaN(); }
-                        }
+                        if (this->arr[le] != F{0}) { x[i] = std::numeric_limits<F>::quiet_NaN(); } // else x[i] remains 0
                     }
                 }
             }
+
             return x;
         }
 
diff --git a/tests/mat_4x4_eigenvalues.cpp b/tests/mat_4x4_eigenvalues.cpp
index 698ca0a..92807bc 100644
--- a/tests/mat_4x4_eigenvalues.cpp
+++ b/tests/mat_4x4_eigenvalues.cpp
@@ -38,6 +38,15 @@ int main()
         }
         std::cout << "\n";
 
+        if ((eigenvalues[0].real() - A[0]) > 1e-10
+            || (eigenvalues[1].real() - A[5]) > 1e-10
+            || (eigenvalues[2].real() - A[10]) > 1e-10
+            || (eigenvalues[3].real() - A[15]) > 1e-10) {
+            std::cout << "Eigenvalues do not seem to be correct\n";
+            --rtn;
+        }
+
+
         // For diagonal matrix, we verify trace and determinant are correct
         double trace = A.trace();
         double det = A.determinant();