diff --git a/AUTHORS b/AUTHORS
index 5ee87d1522..6c5e9274d7 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -233,6 +233,5 @@ BuildTools <anikpatel1322@gmail.com>
 Anik Patel <74193405+Bobingstern@users.noreply.github.com>
 Vrushaket Chaudhari <82214275+vrushaket@users.noreply.github.com>
 Praise Nnamonu <110940850+praisennamonu1@users.noreply.github.com>
-vrushaket <vrushu00@gmail.com>
 
 # Generated by tools/update-authors.js
diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index cdf4172fe0..50edd5cb70 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -1,22 +1,34 @@
 import { factory } from '../../utils/factory.js'
 import { format } from '../../utils/string.js'
 import { createComplexEigs } from './eigs/complexEigs.js'
-import { createRealSymmetric } from './eigs/realSymetric.js'
+import { createRealSymmetric } from './eigs/realSymmetric.js'
 import { typeOf, isNumber, isBigNumber, isComplex, isFraction } from '../../utils/is.js'
 
 const name = 'eigs'
 
 // The absolute state of math.js's dependency system:
-const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolve', 'usolveAll', 'im', 're', 'smaller', 'matrixFromColumns', 'dot']
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, qr, usolve, usolveAll, im, re, smaller, matrixFromColumns, dot }) => {
-  const doRealSymetric = createRealSymmetric({ config, addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplexEigs = createComplexEigs({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, qr, inv, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot })
+const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'size', 'reshape', 'qr', 'usolve', 'usolveAll', 'im', 're', 'smaller', 'matrixFromColumns', 'dot']
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, size, reshape, qr, usolve, usolveAll, im, re, smaller, matrixFromColumns, dot }) => {
+  const doRealSymmetric = createRealSymmetric({ config, addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
+  const doComplexEigs = createComplexEigs({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, size, reshape, qr, inv, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot })
 
   /**
-   * Compute eigenvalues and eigenvectors of a matrix. The eigenvalues are sorted by their absolute value, ascending.
-   * An eigenvalue with multiplicity k will be listed k times. The eigenvectors are returned as columns of a matrix –
-   * the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`).
-   * If the algorithm fails to converge, it will throw an error – in that case, however, you may still find useful information
+   * Compute eigenvalues and eigenvectors of a square matrix.
+   * The eigenvalues are sorted by their absolute value, ascending, and
+   * returned as a vector in the `values` property of the returned project.
+   * An eigenvalue with algebraic multiplicity k will be listed k times, so
+   * that the returned `values` vector always has length equal to the size
+   * of the input matrix.
+   *
+   * The `eigenvectors` property of the return value provides the eigenvectors.
+   * It is an array of plain objects: the `value` property of each gives the
+   * associated eigenvalue, and the `vector` property gives the eigenvector
+   * itself. Note that the same `value` property will occur as many times in
+   * the list provided by `eigenvectors` as the geometric multiplicity of
+   * that value.
+   *
+   * If the algorithm fails to converge, it will throw an error –
+   * in that case, however, you may still find useful information
    * in `err.values` and `err.vectors`.
    *
    * Syntax:
@@ -25,14 +37,15 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
    *
    * Examples:
    *
-   *     const { eigs, multiply, column, transpose } = math
+   *     const { eigs, multiply, column, transpose, matrixFromColumns } = math
    *     const H = [[5, 2.3], [2.3, 1]]
-   *     const ans = eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}
+   *     const ans = eigs(H) // returns {values: [E1,E2...sorted], eigenvectors: [{value: E1, vector: v2}, {value: e, vector: v2}, ...]
    *     const E = ans.values
-   *     const U = ans.vectors
-   *     multiply(H, column(U, 0)) // returns multiply(E[0], column(U, 0))
-   *     const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H
-   *     E[0] == UTxHxU[0][0]  // returns true
+   *     const V = ans.eigenvectors
+   *     multiply(H, V[0].vector)) // returns multiply(E[0], V[0].vector))
+   *     const U = matrixFromColumns(...V.map(obj => obj.vector))
+   *     const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H if possible
+   *     E[0] == UTxHxU[0][0]  // returns true always
    *
    * See also:
    *
@@ -41,62 +54,71 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
    * @param {Array | Matrix} x  Matrix to be diagonalized
    *
    * @param {number | BigNumber} [prec] Precision, default value: 1e-15
-   * @return {{values: Array|Matrix, vectors: Array|Matrix}} Object containing an array of eigenvalues and a matrix with eigenvectors as columns.
+   * @return {{values: Array|Matrix, eigenvectors: Array<EVobj>}} Object containing an array of eigenvalues and an array of {value: number|BigNumber, vector: Array|Matrix} objects.
    *
    */
   return typed('eigs', {
 
-    Array: function (x) {
-      const mat = matrix(x)
-      return computeValuesAndVectors(mat)
-    },
-
+    // The conversion to matrix in the first two implementations,
+    // just to convert back to an array right away in
+    // computeValuesAndVectors, is unfortunate, and should perhaps be
+    // streamlined. It is done because the Matrix object carries some
+    // type information about its entries, and so constructing the matrix
+    // is a roundabout way of doing type detection.
+    Array: function (x) { return doEigs(matrix(x)) },
     'Array, number|BigNumber': function (x, prec) {
-      const mat = matrix(x)
-      return computeValuesAndVectors(mat, prec)
+      return doEigs(matrix(x), prec)
     },
-
     Matrix: function (mat) {
-      const { values, vectors } = computeValuesAndVectors(mat)
-      return {
-        values: matrix(values),
-        vectors: matrix(vectors)
-      }
+      return doEigs(mat, undefined, true)
     },
-
     'Matrix, number|BigNumber': function (mat, prec) {
-      const { values, vectors } = computeValuesAndVectors(mat, prec)
-      return {
-        values: matrix(values),
-        vectors: matrix(vectors)
-      }
+      return doEigs(mat, prec, true)
     }
   })
 
+  function doEigs (mat, prec, matricize = false) {
+    const result = computeValuesAndVectors(mat, prec)
+    if (matricize) {
+      result.values = matrix(result.values)
+      result.eigenvectors = result.eigenvectors.map(({ value, vector }) =>
+        ({ value, vector: matrix(vector) }))
+    }
+    Object.defineProperty(result, 'vectors', {
+      enumerable: false, // to make sure that the eigenvectors can still be
+      // converted to string.
+      get: () => {
+        throw new Error('eigs(M).vectors replaced with eigs(M).eigenvectors')
+      }
+    })
+    return result
+  }
+
   function computeValuesAndVectors (mat, prec) {
     if (prec === undefined) {
       prec = config.epsilon
     }
 
-    const size = mat.size()
+    const arr = mat.toArray() // NOTE: arr is guaranteed to be unaliased
+    // and so safe to modify in place
+    const asize = mat.size()
 
-    if (size.length !== 2 || size[0] !== size[1]) {
-      throw new RangeError('Matrix must be square (size: ' + format(size) + ')')
+    if (asize.length !== 2 || asize[0] !== asize[1]) {
+      throw new RangeError(`Matrix must be square (size: ${format(asize)})`)
     }
 
-    const arr = mat.toArray()
-    const N = size[0]
+    const N = asize[0]
 
     if (isReal(arr, N, prec)) {
-      coerceReal(arr, N)
+      coerceReal(arr, N) // modifies arr by side effect
 
       if (isSymmetric(arr, N, prec)) {
-        const type = coerceTypes(mat, arr, N)
-        return doRealSymetric(arr, N, prec, type)
+        const type = coerceTypes(mat, arr, N) // modifies arr by side effect
+        return doRealSymmetric(arr, N, prec, type)
       }
     }
 
-    const type = coerceTypes(mat, arr, N)
+    const type = coerceTypes(mat, arr, N) // modifies arr by side effect
     return doComplexEigs(arr, N, prec, type)
   }
 
diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js
index 2e337fc207..fc376577f1 100644
--- a/src/function/matrix/eigs/complexEigs.js
+++ b/src/function/matrix/eigs/complexEigs.js
@@ -1,6 +1,6 @@
 import { clone } from '../../../utils/object.js'
 
-export function createComplexEigs ({ addScalar, subtract, flatten, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, inv, qr, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot }) {
+export function createComplexEigs ({ addScalar, subtract, flatten, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, size, reshape, inv, qr, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot }) {
   /**
    * @param {number[][]} arr the matrix to find eigenvalues of
    * @param {number} N size of the matrix
@@ -23,9 +23,9 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     const R = balance(arr, N, prec, type, findVectors)
 
     // R is the row transformation matrix
-    // arr = A' = R A R⁻¹, A is the original matrix
+    // arr = A' = R A R^-1, A is the original matrix
     // (if findVectors is false, R is undefined)
-    // (And so to return to original matrix: A = R⁻¹ arr R)
+    // (And so to return to original matrix: A = R^-1 arr R)
 
     // TODO if magnitudes of elements vary over many orders,
     // move greatest elements to the top left corner
@@ -35,7 +35,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     // updates the transformation matrix R with new row operationsq
     // MODIFIES arr by side effect!
     reduceToHessenberg(arr, N, prec, type, findVectors, R)
-    // still true that original A = R⁻¹ arr R)
+    // still true that original A = R^-1 arr R)
 
     // find eigenvalues
     const { values, C } = iterateUntilTriangular(arr, N, prec, type, findVectors)
@@ -43,17 +43,16 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     // values is the list of eigenvalues, C is the column
     // transformation matrix that transforms arr, the hessenberg
     // matrix, to upper triangular
-    // (So U = C⁻¹ arr C and the relationship between current arr
+    // (So U = C^-1 arr C and the relationship between current arr
     // and original A is unchanged.)
 
-    let vectors
+    let eigenvectors
 
     if (findVectors) {
-      vectors = findEigenvectors(arr, N, C, R, values, prec, type)
-      vectors = matrixFromColumns(...vectors)
+      eigenvectors = findEigenvectors(arr, N, C, R, values, prec, type)
     }
 
-    return { values, vectors }
+    return { values, eigenvectors }
   }
 
   /**
@@ -255,7 +254,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
 
     // The Francis Algorithm
     // The core idea of this algorithm is that doing successive
-    // A' = Q⁺AQ transformations will eventually converge to block-
+    // A' = QtAQ transformations will eventually converge to block-
     // upper-triangular with diagonal blocks either 1x1 or 2x2.
     // The Q here is the one from the QR decomposition, A = QR.
     // Since the eigenvalues of a block-upper-triangular matrix are
@@ -277,7 +276,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     // N×N matrix describing the overall transformation done during the QR algorithm
     let Qtotal = findVectors ? diag(Array(N).fill(one)) : undefined
 
-    // n×n matrix describing the QR transformations done since last convergence
+    // nxn matrix describing the QR transformations done since last convergence
     let Qpartial = findVectors ? diag(Array(n).fill(one)) : undefined
 
     // last eigenvalue converged before this many steps
@@ -290,7 +289,12 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
 
       // Perform the factorization
 
-      const k = 0 // TODO set close to an eigenvalue
+      const k = arr[n - 1][n - 1] // TODO this is apparently a somewhat
+      // old-fashioned choice; ideally set close to an eigenvalue, or
+      // perhaps better yet switch to the implicit QR version that is sometimes
+      // specifically called the "Francis algorithm" that is alluded to
+      // in the following TODO. (Or perhaps we switch to an independently
+      // optimized third-party package for the linear algebra operations...)
 
       for (let i = 0; i < n; i++) {
         arr[i][i] = subtract(arr[i][i], k)
@@ -350,7 +354,6 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
           ))
           inflateMatrix(Qpartial, N)
           Qtotal = multiply(Qtotal, Qpartial)
-
           if (n > 2) {
             Qpartial = diag(Array(n - 2).fill(one))
           }
@@ -413,18 +416,18 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     const uniqueValues = []
     const multiplicities = []
 
-    for (const λ of values) {
-      const i = indexOf(uniqueValues, λ, equal)
+    for (const lambda of values) {
+      const i = indexOf(uniqueValues, lambda, equal)
 
       if (i === -1) {
-        uniqueValues.push(λ)
+        uniqueValues.push(lambda)
         multiplicities.push(1)
       } else {
         multiplicities[i] += 1
       }
     }
 
-    // find eigenvectors by solving U − λE = 0
+    // find eigenvectors by solving U − lambdaE = 0
     // TODO replace with an iterative eigenvector algorithm
     // (this one might fail for imprecise eigenvalues)
 
@@ -433,26 +436,19 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     const b = Array(N).fill(zero)
     const E = diag(Array(N).fill(one))
 
-    // eigenvalues for which usolve failed (due to numerical error)
-    const failedLambdas = []
-
     for (let i = 0; i < len; i++) {
-      const λ = uniqueValues[i]
-      const S = subtract(U, multiply(λ, E)) // the characteristic matrix
+      const lambda = uniqueValues[i]
+      const S = subtract(U, multiply(lambda, E)) // the characteristic matrix
 
       let solutions = usolveAll(S, b)
       solutions.shift() // ignore the null vector
 
       // looks like we missed something, try inverse iteration
+      // But if that fails, just presume that the original matrix truly
+      // was defective.
       while (solutions.length < multiplicities[i]) {
         const approxVec = inverseIterate(S, N, solutions, prec, type)
-
-        if (approxVec == null) {
-          // no more vectors were found
-          failedLambdas.push(λ)
-          break
-        }
-
+        if (approxVec === null) { break } // no more vectors were found
         solutions.push(approxVec)
       }
 
@@ -460,14 +456,8 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
       const correction = multiply(inv(R), C)
       solutions = solutions.map(v => multiply(correction, v))
 
-      vectors.push(...solutions.map(v => flatten(v)))
-    }
-
-    if (failedLambdas.length !== 0) {
-      const err = new Error('Failed to find eigenvectors for the following eigenvalues: ' + failedLambdas.join(', '))
-      err.values = values
-      err.vectors = vectors
-      throw err
+      vectors.push(
+        ...solutions.map(v => ({ value: lambda, vector: flatten(v) })))
     }
 
     return vectors
@@ -478,7 +468,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
    * @return {[number,number]}
    */
   function eigenvalues2x2 (a, b, c, d) {
-    // λ± = ½ trA ± ½ √( tr²A - 4 detA )
+    // lambda_+- = 1/2 trA +- 1/2 sqrt( tr^2 A - 4 detA )
     const trA = addScalar(a, d)
     const detA = subtract(multiplyScalar(a, d), multiplyScalar(b, c))
     const x = multiplyScalar(trA, 0.5)
@@ -489,7 +479,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
 
   /**
    * For an 2x2 matrix compute the transformation matrix S,
-   * so that SAS⁻¹ is an upper triangular matrix
+   * so that SAS^-1 is an upper triangular matrix
    * @return {[[number,number],[number,number]]}
    * @see https://math.berkeley.edu/~ogus/old/Math_54-05/webfoils/jordan.pdf
    * @see http://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html
@@ -514,24 +504,22 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     }
 
     // matrix is not diagonalizable
-    // compute off-diagonal elements of N = A - λI
-    // N₁₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ )
-    // N₁₂ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ )
-
+    // compute diagonal elements of N = A - lambdaI
     const na = subtract(a, l1)
-    const nb = subtract(b, l1)
-    const nc = subtract(c, l1)
     const nd = subtract(d, l1)
 
-    if (smaller(abs(nb), prec)) {
-      return [[na, one], [nc, zero]]
+    // col(N,2) = 0  implies  S = ( col(N,1), e_1 )
+    // col(N,2) != 0 implies  S = ( col(N,2), e_2 )
+
+    if (smaller(abs(b), prec) && smaller(abs(nd), prec)) {
+      return [[na, one], [c, zero]]
     } else {
-      return [[nb, zero], [nd, one]]
+      return [[b, zero], [nd, one]]
     }
   }
 
   /**
-   * Enlarge the matrix from n×n to N×N, setting the new
+   * Enlarge the matrix from nxn to NxN, setting the new
    * elements to 1 on diagonal and 0 elsewhere
    */
   function inflateMatrix (arr, N) {
@@ -614,12 +602,19 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
 
     // you better choose a random vector before I count to five
     let i = 0
-    while (true) {
+    for (; i < 5; ++i) {
       b = randomOrthogonalVector(N, orthog, type)
-      b = usolve(A, b)
-
+      try {
+        b = usolve(A, b)
+      } catch {
+        // That direction didn't work, likely because the original matrix
+        // was defective. But still make the full number of tries...
+        continue
+      }
       if (larger(norm(b), largeNum)) { break }
-      if (++i >= 5) { return null }
+    }
+    if (i >= 5) {
+      return null // couldn't find any orthogonal vector in the image
     }
 
     // you better converge before I count to ten
@@ -664,8 +659,10 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
    * Project vector v to the orthogonal complement of an array of vectors
    */
   function orthogonalComplement (v, orthog) {
-    for (const w of orthog) {
-      // v := v − (w, v)/∥w∥² w
+    const vectorShape = size(v)
+    for (let w of orthog) {
+      w = reshape(w, vectorShape) // make sure this is just a vector computation
+      // v := v − (w, v)/|w|^2 w
       v = subtract(v, multiply(divideScalar(dot(w, v), dot(w, w)), w))
     }
 
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymmetric.js
similarity index 86%
rename from src/function/matrix/eigs/realSymetric.js
rename to src/function/matrix/eigs/realSymmetric.js
index 57c215f2c6..7e24d406da 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymmetric.js
@@ -27,7 +27,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     let Sij = new Array(N)
     // Sij is Identity Matrix
     for (let i = 0; i < N; i++) {
-      Sij[i] = createArray(N, 0)
+      Sij[i] = Array(N).fill(0)
       Sij[i][i] = 1.0
     }
     // initial error
@@ -40,7 +40,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       Sij = Sij1(Sij, psi, i, j)
       Vab = getAij(x)
     }
-    const Ei = createArray(N, 0) // eigenvalues
+    const Ei = Array(N).fill(0) // eigenvalues
     for (let i = 0; i < N; i++) {
       Ei[i] = x[i][i]
     }
@@ -55,7 +55,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     let Sij = new Array(N)
     // Sij is Identity Matrix
     for (let i = 0; i < N; i++) {
-      Sij[i] = createArray(N, 0)
+      Sij[i] = Array(N).fill(0)
       Sij[i][i] = 1.0
     }
     // initial error
@@ -68,7 +68,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       Sij = Sij1Big(Sij, psi, i, j)
       Vab = getAijBig(x)
     }
-    const Ei = createArray(N, 0) // eigenvalues
+    const Ei = Array(N).fill(0) // eigenvalues
     for (let i = 0; i < N; i++) {
       Ei[i] = x[i][i]
     }
@@ -101,8 +101,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const N = Sij.length
     const c = Math.cos(theta)
     const s = Math.sin(theta)
-    const Ski = createArray(N, 0)
-    const Skj = createArray(N, 0)
+    const Ski = Array(N).fill(0)
+    const Skj = Array(N).fill(0)
     for (let k = 0; k < N; k++) {
       Ski[k] = c * Sij[k][i] - s * Sij[k][j]
       Skj[k] = s * Sij[k][i] + c * Sij[k][j]
@@ -118,8 +118,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const N = Sij.length
     const c = cos(theta)
     const s = sin(theta)
-    const Ski = createArray(N, bignumber(0))
-    const Skj = createArray(N, bignumber(0))
+    const Ski = Array(N).fill(bignumber(0))
+    const Skj = Array(N).fill(bignumber(0))
     for (let k = 0; k < N; k++) {
       Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]))
       Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]))
@@ -138,8 +138,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const s = bignumber(sin(theta))
     const c2 = multiplyScalar(c, c)
     const s2 = multiplyScalar(s, s)
-    const Aki = createArray(N, bignumber(0))
-    const Akj = createArray(N, bignumber(0))
+    const Aki = Array(N).fill(bignumber(0))
+    const Akj = Array(N).fill(bignumber(0))
     // 2cs Hij
     const csHij = multiply(bignumber(2), c, s, Hij[i][j])
     //  Aii
@@ -174,8 +174,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const s = Math.sin(theta)
     const c2 = c * c
     const s2 = s * s
-    const Aki = createArray(N, 0)
-    const Akj = createArray(N, 0)
+    const Aki = Array(N).fill(0)
+    const Akj = Array(N).fill(0)
     //  Aii
     const Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j]
     const Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j]
@@ -237,10 +237,10 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
   function sorting (E, S) {
     const N = E.length
     const values = Array(N)
-    const vectors = Array(N)
+    const vecs = Array(N)
 
     for (let k = 0; k < N; k++) {
-      vectors[k] = Array(N)
+      vecs[k] = Array(N)
     }
     for (let i = 0; i < N; i++) {
       let minID = 0
@@ -253,29 +253,12 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       }
       values[i] = E.splice(minID, 1)[0]
       for (let k = 0; k < N; k++) {
-        vectors[k][i] = S[k][minID]
+        vecs[i][k] = S[k][minID]
         S[k].splice(minID, 1)
       }
     }
-
-    return { values, vectors }
-  }
-
-  /**
-   * Create an array of a certain size and fill all items with an initial value
-   * @param {number} size
-   * @param {number} value
-   * @return {number[]}
-   */
-  function createArray (size, value) {
-    // TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it)
-    const array = new Array(size)
-
-    for (let i = 0; i < size; i++) {
-      array[i] = value
-    }
-
-    return array
+    const eigenvectors = vecs.map((vector, i) => ({ value: values[i], vector }))
+    return { values, eigenvectors }
   }
 
   return main
diff --git a/src/utils/number.js b/src/utils/number.js
index 730f2cf7ac..fd72b57c04 100644
--- a/src/utils/number.js
+++ b/src/utils/number.js
@@ -623,7 +623,7 @@ export function nearlyEqual (x, y, epsilon) {
   if (isFinite(x) && isFinite(y)) {
     // check numbers are very close, needed when comparing numbers near zero
     const diff = Math.abs(x - y)
-    if (diff < DBL_EPSILON) {
+    if (diff <= DBL_EPSILON) {
       return true
     } else {
       // use relative error
diff --git a/test/typescript-tests/testTypes.ts b/test/typescript-tests/testTypes.ts
index 1a777c71c4..23a1a34fd9 100644
--- a/test/typescript-tests/testTypes.ts
+++ b/test/typescript-tests/testTypes.ts
@@ -1282,6 +1282,17 @@ Matrices examples
     assert.strictEqual(math.matrix([1, 2, 3]) instanceof math.Matrix, true)
   }
 
+  // Eigenvalues and eigenvectors
+  {
+    const D = [
+      [1, 1],
+      [0, 1]
+    ]
+    const eig = math.eigs(D)
+    assert.ok(math.deepEqual(eig.values, [1, 1]))
+    assert.deepStrictEqual(eig.eigenvectors, [{ value: 1, vector: [1, 0] }])
+  }
+
   // Fourier transform and inverse
   {
     assert.ok(
diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js
index f907b1bd61..f77862d651 100644
--- a/test/unit-tests/function/matrix/eigs.test.js
+++ b/test/unit-tests/function/matrix/eigs.test.js
@@ -1,9 +1,14 @@
 import assert from 'assert'
 import math from '../../../../src/defaultInstance.js'
 import approx from '../../../../tools/approx.js'
-const { eigs, add, complex, divide, exp, fraction, matrix, multiply, size, transpose, bignumber: bignum, zeros, Matrix, Complex } = math
+const { eigs, add, complex, divide, exp, fraction, matrix, matrixFromColumns, multiply, abs, size, transpose, bignumber: bignum, zeros, Matrix, Complex } = math
 
 describe('eigs', function () {
+  // helper to examine eigenvectors
+  function testEigenvectors (soln, predicate) {
+    soln.eigenvectors.forEach((ev, i) => predicate(ev.vector, i))
+  }
+
   it('only accepts a square matrix', function () {
     assert.throws(function () { eigs(matrix([[1, 2, 3], [4, 5, 6]])) }, /Matrix must be square/)
     assert.throws(function () { eigs([[1, 2, 3], [4, 5, 6]]) }, /Matrix must be square/)
@@ -16,23 +21,30 @@ describe('eigs', function () {
   it('follows aiao-mimo', function () {
     const realSymArray = eigs([[1, 0], [0, 1]])
     assert(Array.isArray(realSymArray.values) && typeof realSymArray.values[0] === 'number')
-    assert(Array.isArray(realSymArray.vectors) && typeof realSymArray.vectors[0][0] === 'number')
+    testEigenvectors(realSymArray, vector => assert(Array.isArray(vector)))
+    assert(typeof realSymArray.eigenvectors[0].vector[0] === 'number')
 
     const genericArray = eigs([[0, 1], [-1, 0]])
     assert(Array.isArray(genericArray.values) && genericArray.values[0] instanceof Complex)
-    assert(Array.isArray(genericArray.vectors) && genericArray.vectors[0][0] instanceof Complex)
+    testEigenvectors(genericArray,
+      vector => assert(Array.isArray(vector) && vector[0] instanceof Complex)
+    )
 
     const realSymMatrix = eigs(matrix([[1, 0], [0, 1]]))
     assert(realSymMatrix.values instanceof Matrix)
     assert.deepStrictEqual(size(realSymMatrix.values), matrix([2]))
-    assert(realSymMatrix.vectors instanceof Matrix)
-    assert.deepStrictEqual(size(realSymMatrix.vectors), matrix([2, 2]))
+    testEigenvectors(realSymMatrix, vector => {
+      assert(vector instanceof Matrix)
+      assert.deepStrictEqual(size(vector), matrix([2]))
+    })
 
     const genericMatrix = eigs(matrix([[0, 1], [-1, 0]]))
     assert(genericMatrix.values instanceof Matrix)
     assert.deepStrictEqual(size(genericMatrix.values), matrix([2]))
-    assert(genericMatrix.vectors instanceof Matrix)
-    assert.deepStrictEqual(size(genericMatrix.vectors), matrix([2, 2]))
+    testEigenvectors(genericMatrix, vector => {
+      assert(vector instanceof Matrix)
+      assert.deepStrictEqual(size(vector), matrix([2]))
+    })
   })
 
   it('only accepts a matrix with valid element type', function () {
@@ -117,8 +129,7 @@ describe('eigs', function () {
 
     // inverse iteration is stochastic, check it multiple times
     for (let i = 0; i < 5; i++) {
-      const { vectors } = eigs(m)
-      const eigenRows = transpose(vectors)
+      const eigenRows = eigs(m).eigenvectors.map(obj => obj.vector)
       // if we scale each row to the expected scale, they should match
       for (let j = 0; j < 5; j++) {
         approx.deepEqual(divide(eigenRows[i], eigenRows[i][oneIndex[i]]),
@@ -136,7 +147,11 @@ describe('eigs', function () {
       [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]]
     const ans = eigs(H)
     const E = ans.values
-    const V = ans.vectors
+    testEigenvectors(ans,
+      (v, j) => approx.deepEqual(multiply(E[j], v), multiply(H, v))
+    )
+    const Vcols = ans.eigenvectors.map(obj => obj.vector)
+    const V = matrixFromColumns(...Vcols)
     const VtHV = multiply(transpose(V), H, V)
     const Ei = Array(H.length)
     for (let i = 0; i < H.length; i++) {
@@ -151,12 +166,10 @@ describe('eigs', function () {
     const cnt = 0.1
     const Ath = multiply(exp(multiply(complex(0, 1), -cnt)), A)
     const Hth = divide(add(Ath, transpose(Ath)), 2)
-    const { values, vectors } = eigs(Hth)
-    const R = transpose(vectors) // rows are eigenvectors
-    for (const i of [0, 1, 2]) {
-      const v = R[i]
-      approx.deepEqual(multiply(Hth, v), multiply(values[i], v))
-    }
+    const example = eigs(Hth)
+    testEigenvectors(example, (v, i) =>
+      approx.deepEqual(multiply(Hth, v), multiply(example.values[i], v))
+    )
   })
 
   it('supports fractions', function () {
@@ -169,6 +182,61 @@ describe('eigs', function () {
     )
   })
 
+  it('handles some 2x2 defective matrices', function () {
+    const check = eigs([[2.0, 1.0], [0.0, 2.0]]) // Test case from #2879
+    assert.deepStrictEqual(check, {
+      values: [2, 2],
+      eigenvectors: [{ value: 2, vector: [1, 0] }]
+    })
+    const fromWeb = eigs([[-2, 1], [-1, 0]]) // https://ocw.mit.edu/courses/18-03sc-differential-equations-fall-2011/051316d5fa93f560934d3e410f8d153d_MIT18_03SCF11_s33_8text.pdf
+    assert.strictEqual(fromWeb.eigenvectors.length, 1)
+    const vec = fromWeb.eigenvectors[0].vector
+    approx.equal(vec[0], vec[1])
+  })
+
+  it('handles a 3x3 defective matrix', function () {
+    const fromWeb = eigs([[2, -5, 0], [0, 2, 0], [-1, 4, 1]]) // https://math.libretexts.org/Bookshelves/Differential_Equations/Differential_Equations_for_Engineers_(Lebl)/3%3A_Systems_of_ODEs/3.7%3A_Multiple_Eigenvalues
+    assert.strictEqual(fromWeb.eigenvectors.length, 2)
+    const ev = fromWeb.eigenvectors
+    approx.equal(ev[0].value, 1)
+    approx.equal(ev[1].value, 2)
+    approx.equal(ev[0].vector[0], 0)
+    approx.equal(ev[0].vector[1], 0)
+    assert.ok(abs(ev[0].vector[2]) > math.config.epsilon)
+    approx.equal(ev[1].vector[0], -ev[1].vector[2])
+    approx.equal(ev[1].vector[1], 0)
+    const web2 = eigs([[1, 1, 0], [0, 1, 2], [0, 0, 3]]) // https://www2.math.upenn.edu/~moose/240S2013/slides7-31.pdf
+    assert.strictEqual(web2.eigenvectors.length, 2)
+    const ev2 = web2.eigenvectors
+    assert.strictEqual(ev2[0].value, 1)
+    assert.strictEqual(ev2[1].value, 3)
+    assert.strictEqual(ev2[0].vector[1], 0)
+    assert.strictEqual(ev2[0].vector[2], 0)
+    assert.ok(abs(ev2[0].vector[0]) > math.config.epsilon)
+    assert.strictEqual(ev2[1].vector[1], ev2[1].vector[2])
+    approx.equal(ev2[1].vector[1], 2 * ev2[1].vector[0])
+  })
+
+  it('accepts a precision argument', function () {
+    // The following is a matrix with an algebraically triple eigenvalue
+    // equal to 2 which has a unique eigenvector (up to scale, of course).
+    // It is from https://web.uvic.ca/~tbazett/diffyqs/sec_multeigen.html
+    // The iterative eigenvalue calculation currently being used has a
+    // great deal of difficulty converging. We can use a fine precision,
+    // but it still doesn't produce good eigenvalues. Hopefully someday
+    // we'll be able to get closer.
+    const difficult = [[2, 0, 0], [-1, -1, 9], [0, -1, 5]]
+    const poor = eigs(difficult, 1e-14)
+    assert.strictEqual(poor.values.length, 3)
+    approx.deepEqual(poor.values, [2, 2, 2], 6e-6)
+    // Note the eigenvectors are junk, so we don't test them. The function
+    // eigs thinks there are three of them, for example. Hopefully some
+    // future iteration of mathjs will be able to discover there is really
+    // only one.
+    const poorm = eigs(matrix(difficult), 1e-14)
+    assert.deepStrictEqual(poorm.values.size(), [3])
+  })
+
   it('diagonalizes matrix with bigNumber', function () {
     const x = [[bignum(1), bignum(0)], [bignum(0), bignum(1)]]
     approx.deepEqual(eigs(x).values, [bignum(1), bignum(1)])
@@ -184,7 +252,8 @@ describe('eigs', function () {
       [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]])
     const ans = eigs(H)
     const E = ans.values
-    const V = ans.vectors
+    const Vcols = ans.eigenvectors.map(obj => obj.vector)
+    const V = matrixFromColumns(...Vcols)
     const VtHV = multiply(transpose(V), H, V)
     const Ei = Array(H.length)
     for (let i = 0; i < H.length; i++) {
diff --git a/types/index.d.ts b/types/index.d.ts
index c14a05d3f9..a013d1c2d5 100644
--- a/types/index.d.ts
+++ b/types/index.d.ts
@@ -1753,9 +1753,8 @@ declare namespace math {
      * Compute eigenvalues and eigenvectors of a matrix.
      * The eigenvalues are sorted by their absolute value, ascending.
      * An eigenvalue with multiplicity k will be listed k times.
-     * The eigenvectors are returned as columns of a matrix – the eigenvector
-     * that belongs to the j-th eigenvalue in the list (eg. values[j]) is the
-     * j-th column (eg. column(vectors, j)). If the algorithm fails to converge,
+     * The eigenvectors are returned as an array of objects, each with a
+     * `value` and a `vector`. If the algorithm fails to converge,
      * it will throw an error – in that case, however, you may still find useful
      * information in err.values and err.vectors
      * @param x Matrix to be diagonalized
@@ -1765,7 +1764,13 @@ declare namespace math {
     eigs(
       x: MathCollection,
       prec?: number | BigNumber
-    ): { values: MathCollection; vectors: MathCollection }
+    ): {
+      values: MathCollection
+      eigenvectors: {
+        value: number | BigNumber
+        vector: MathCollection
+      }[]
+    }
 
     /**
      * Compute the matrix exponential, expm(A) = e^A. The matrix must be