diff --git a/src/expression/embeddedDocs/embeddedDocs.js b/src/expression/embeddedDocs/embeddedDocs.js
index 5c8bfb2267..840555c970 100644
--- a/src/expression/embeddedDocs/embeddedDocs.js
+++ b/src/expression/embeddedDocs/embeddedDocs.js
@@ -167,11 +167,13 @@ import { addDocs } from './function/arithmetic/add'
 import { absDocs } from './function/arithmetic/abs'
 import { qrDocs } from './function/algebra/qr'
 import { usolveDocs } from './function/algebra/usolve'
+import { usolveAllDocs } from './function/algebra/usolveAll'
 import { sluDocs } from './function/algebra/slu'
 import { rationalizeDocs } from './function/algebra/rationalize'
 import { simplifyDocs } from './function/algebra/simplify'
 import { lupDocs } from './function/algebra/lup'
 import { lsolveDocs } from './function/algebra/lsolve'
+import { lsolveAllDocs } from './function/algebra/lsolveAll'
 import { derivativeDocs } from './function/algebra/derivative'
 import { versionDocs } from './constants/version'
 import { trueDocs } from './constants/true'
@@ -310,12 +312,14 @@ export const embeddedDocs = {
   // functions - algebra
   derivative: derivativeDocs,
   lsolve: lsolveDocs,
+  lsolveAll: lsolveAllDocs,
   lup: lupDocs,
   lusolve: lusolveDocs,
   simplify: simplifyDocs,
   rationalize: rationalizeDocs,
   slu: sluDocs,
   usolve: usolveDocs,
+  usolveAll: usolveAllDocs,
   qr: qrDocs,
 
   // functions - arithmetic
diff --git a/src/expression/embeddedDocs/function/algebra/lsolve.js b/src/expression/embeddedDocs/function/algebra/lsolve.js
index 0cf27bc9d6..7fc0474e6b 100644
--- a/src/expression/embeddedDocs/function/algebra/lsolve.js
+++ b/src/expression/embeddedDocs/function/algebra/lsolve.js
@@ -5,13 +5,13 @@ export const lsolveDocs = {
     'x=lsolve(L, b)'
   ],
   description:
-  'Solves the linear system L * x = b where L is an [n x n] lower triangular matrix and b is a [n] column vector.',
+  'Finds one solution of the linear system L * x = b where L is an [n x n] lower triangular matrix and b is a [n] column vector.',
   examples: [
     'a = [-2, 3; 2, 1]',
     'b = [11, 9]',
     'x = lsolve(a, b)'
   ],
   seealso: [
-    'lup', 'lusolve', 'usolve', 'matrix', 'sparse'
+    'lsolveAll', 'lup', 'lusolve', 'usolve', 'matrix', 'sparse'
   ]
 }
diff --git a/src/expression/embeddedDocs/function/algebra/lsolveAll.js b/src/expression/embeddedDocs/function/algebra/lsolveAll.js
new file mode 100644
index 0000000000..e5ff44cf24
--- /dev/null
+++ b/src/expression/embeddedDocs/function/algebra/lsolveAll.js
@@ -0,0 +1,17 @@
+export const lsolveAllDocs = {
+  name: 'lsolveAll',
+  category: 'Algebra',
+  syntax: [
+    'x=lsolveAll(L, b)'
+  ],
+  description:
+  'Finds all solutions of the linear system L * x = b where L is an [n x n] lower triangular matrix and b is a [n] column vector.',
+  examples: [
+    'a = [-2, 3; 2, 1]',
+    'b = [11, 9]',
+    'x = lsolve(a, b)'
+  ],
+  seealso: [
+    'lsolve', 'lup', 'lusolve', 'usolve', 'matrix', 'sparse'
+  ]
+}
diff --git a/src/expression/embeddedDocs/function/algebra/usolve.js b/src/expression/embeddedDocs/function/algebra/usolve.js
index 7c0903a1de..9a22974a37 100644
--- a/src/expression/embeddedDocs/function/algebra/usolve.js
+++ b/src/expression/embeddedDocs/function/algebra/usolve.js
@@ -5,11 +5,11 @@ export const usolveDocs = {
     'x=usolve(U, b)'
   ],
   description:
-  'Solves the linear system U * x = b where U is an [n x n] upper triangular matrix and b is a [n] column vector.',
+  'Finds one solution of the linear system U * x = b where U is an [n x n] upper triangular matrix and b is a [n] column vector.',
   examples: [
     'x=usolve(sparse([1, 1, 1, 1; 0, 1, 1, 1; 0, 0, 1, 1; 0, 0, 0, 1]), [1; 2; 3; 4])'
   ],
   seealso: [
-    'lup', 'lusolve', 'lsolve', 'matrix', 'sparse'
+    'usolveAll', 'lup', 'lusolve', 'lsolve', 'matrix', 'sparse'
   ]
 }
diff --git a/src/expression/embeddedDocs/function/algebra/usolveAll.js b/src/expression/embeddedDocs/function/algebra/usolveAll.js
new file mode 100644
index 0000000000..c2131c1be1
--- /dev/null
+++ b/src/expression/embeddedDocs/function/algebra/usolveAll.js
@@ -0,0 +1,15 @@
+export const usolveAllDocs = {
+  name: 'usolveAll',
+  category: 'Algebra',
+  syntax: [
+    'x=usolve(U, b)'
+  ],
+  description:
+  'Finds all solutions of the linear system U * x = b where U is an [n x n] upper triangular matrix and b is a [n] column vector.',
+  examples: [
+    'x=usolve(sparse([1, 1, 1, 1; 0, 1, 1, 1; 0, 0, 1, 1; 0, 0, 0, 1]), [1; 2; 3; 4])'
+  ],
+  seealso: [
+    'usolve', 'lup', 'lusolve', 'lsolve', 'matrix', 'sparse'
+  ]
+}
diff --git a/src/factoriesAny.js b/src/factoriesAny.js
index de60e761aa..a87428fd74 100644
--- a/src/factoriesAny.js
+++ b/src/factoriesAny.js
@@ -103,6 +103,8 @@ export { createDotPow } from './function/arithmetic/dotPow'
 export { createDotDivide } from './function/arithmetic/dotDivide'
 export { createLsolve } from './function/algebra/solver/lsolve'
 export { createUsolve } from './function/algebra/solver/usolve'
+export { createLsolveAll } from './function/algebra/solver/lsolveAll'
+export { createUsolveAll } from './function/algebra/solver/usolveAll'
 export { createLeftShift } from './function/bitwise/leftShift'
 export { createRightArithShift } from './function/bitwise/rightArithShift'
 export { createRightLogShift } from './function/bitwise/rightLogShift'
diff --git a/src/function/algebra/solver/lsolve.js b/src/function/algebra/solver/lsolve.js
index d1e5b024dc..eb1623930f 100644
--- a/src/function/algebra/solver/lsolve.js
+++ b/src/function/algebra/solver/lsolve.js
@@ -16,7 +16,7 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
   const solveValidation = createSolveValidation({ DenseMatrix })
 
   /**
-   * Solves the linear equation system by forwards substitution. Matrix must be a lower triangular matrix.
+   * Finds one solution of a linear equation system by forwards substitution. Matrix must be a lower triangular matrix. Throws an error if there's no solution.
    *
    * `L * x = b`
    *
@@ -32,7 +32,7 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *
    * See also:
    *
-   *    lup, slu, usolve, lusolve
+   *    lsolveAll, lup, slu, usolve, lusolve
    *
    * @param {Matrix, Array} L       A N x N matrix or array (L)
    * @param {Matrix, Array} b       A column vector with the b values
@@ -42,21 +42,16 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
   return typed(name, {
 
     'SparseMatrix, Array | Matrix': function (m, b) {
-      // process matrix
       return _sparseForwardSubstitution(m, b)
     },
 
     'DenseMatrix, Array | Matrix': function (m, b) {
-      // process matrix
       return _denseForwardSubstitution(m, b)
     },
 
     'Array, Array | Matrix': function (a, b) {
-      // create dense matrix from array
       const m = matrix(a)
-      // use matrix implementation
       const r = _denseForwardSubstitution(m, b)
-      // result
       return r.valueOf()
     }
   })
@@ -64,45 +59,44 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
   function _denseForwardSubstitution (m, b) {
     // validate matrix and vector, return copy of column vector b
     b = solveValidation(m, b, true)
-    // column vector data
     const bdata = b._data
-    // rows & columns
+
     const rows = m._size[0]
     const columns = m._size[1]
+
     // result
     const x = []
-    // data
-    const data = m._data
-    // forward solve m * x = b, loop columns
+
+    const mdata = m._data
+
+    // loop columns
     for (let j = 0; j < columns; j++) {
-      // b[j]
       const bj = bdata[j][0] || 0
-      // x[j]
       let xj
-      // forward substitution (outer product) avoids inner looping when bj === 0
+
       if (!equalScalar(bj, 0)) {
-        // value @ [j, j]
-        const vjj = data[j][j]
-        // check vjj
+        // non-degenerate row, find solution
+
+        const vjj = mdata[j][j]
+
         if (equalScalar(vjj, 0)) {
-          // system cannot be solved
           throw new Error('Linear system cannot be solved since matrix is singular')
         }
-        // calculate xj
+
         xj = divideScalar(bj, vjj)
+
         // loop rows
         for (let i = j + 1; i < rows; i++) {
-          // update copy of b
-          bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, data[i][j]))]
+          bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))]
         }
       } else {
-        // zero @ j
+        // degenerate row, we can choose any value
         xj = 0
       }
-      // update x
+
       x[j] = [xj]
     }
-    // return vector
+
     return new DenseMatrix({
       data: x,
       size: [rows, 1]
@@ -112,68 +106,68 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
   function _sparseForwardSubstitution (m, b) {
     // validate matrix and vector, return copy of column vector b
     b = solveValidation(m, b, true)
-    // column vector data
+
     const bdata = b._data
-    // rows & columns
+
     const rows = m._size[0]
     const columns = m._size[1]
-    // matrix arrays
+
     const values = m._values
     const index = m._index
     const ptr = m._ptr
-    // vars
-    let i, k
+
     // result
     const x = []
-    // forward solve m * x = b, loop columns
+
+    // loop columns
     for (let j = 0; j < columns; j++) {
-      // b[j]
       const bj = bdata[j][0] || 0
-      // forward substitution (outer product) avoids inner looping when bj === 0
+
       if (!equalScalar(bj, 0)) {
-        // value @ [j, j]
+        // non-degenerate row, find solution
+
         let vjj = 0
-        // lower triangular matrix values & index (column j)
-        const jvalues = []
-        const jindex = []
-        // last index in column
-        let l = ptr[j + 1]
-        // values in column, find value @ [j, j]
-        for (k = ptr[j]; k < l; k++) {
-          // row
-          i = index[k]
+        // matrix values & indices (column j)
+        const jValues = []
+        const jIndices = []
+
+        // first and last index in the column
+        const firstIndex = ptr[j]
+        const lastIndex = ptr[j + 1]
+
+        // values in column, find value at [j, j]
+        for (let k = firstIndex; k < lastIndex; k++) {
+          const i = index[k]
+
           // check row (rows are not sorted!)
           if (i === j) {
-            // update vjj
             vjj = values[k]
           } else if (i > j) {
             // store lower triangular
-            jvalues.push(values[k])
-            jindex.push(i)
+            jValues.push(values[k])
+            jIndices.push(i)
           }
         }
-        // at this point we must have a value @ [j, j]
+
+        // at this point we must have a value in vjj
         if (equalScalar(vjj, 0)) {
-          // system cannot be solved, there is no value @ [j, j]
           throw new Error('Linear system cannot be solved since matrix is singular')
         }
-        // calculate xj
+
         const xj = divideScalar(bj, vjj)
-        // loop lower triangular
-        for (k = 0, l = jindex.length; k < l; k++) {
-          // row
-          i = jindex[k]
-          // update copy of b
-          bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, jvalues[k]))]
+
+        for (let k = 0, l = jIndices.length; k < l; k++) {
+          const i = jIndices[k]
+          bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, jValues[k]))]
         }
-        // update x
+
         x[j] = [xj]
       } else {
-        // update x
+        // degenerate row, we can choose any value
         x[j] = [0]
       }
     }
-    // return vector
+
     return new DenseMatrix({
       data: x,
       size: [rows, 1]
diff --git a/src/function/algebra/solver/lsolveAll.js b/src/function/algebra/solver/lsolveAll.js
new file mode 100644
index 0000000000..7a4c825785
--- /dev/null
+++ b/src/function/algebra/solver/lsolveAll.js
@@ -0,0 +1,197 @@
+import { factory } from '../../../utils/factory'
+import { createSolveValidation } from './utils/solveValidation'
+
+const name = 'lsolveAll'
+const dependencies = [
+  'typed',
+  'matrix',
+  'divideScalar',
+  'multiplyScalar',
+  'subtract',
+  'equalScalar',
+  'DenseMatrix'
+]
+
+export const createLsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, divideScalar, multiplyScalar, subtract, equalScalar, DenseMatrix }) => {
+  const solveValidation = createSolveValidation({ DenseMatrix })
+
+  /**
+   * Finds all solutions of a linear equation system by forwards substitution. Matrix must be a lower triangular matrix.
+   *
+   * `L * x = b`
+   *
+   * Syntax:
+   *
+   *    math.lsolve(L, b)
+   *
+   * Examples:
+   *
+   *    const a = [[-2, 3], [2, 1]]
+   *    const b = [11, 9]
+   *    const x = lsolve(a, b)  // [ [[-5.5], [20]] ]
+   *
+   * See also:
+   *
+   *    lsolve, lup, slu, usolve, lusolve
+   *
+   * @param {Matrix, Array} L       A N x N matrix or array (L)
+   * @param {Matrix, Array} b       A column vector with the b values
+   *
+   * @return {DenseMatrix[] | Array[]}  An array of affine-independent column vectors (x) that solve the linear system
+   */
+  return typed(name, {
+
+    'SparseMatrix, Array | Matrix': function (m, b) {
+      return _sparseForwardSubstitution(m, b)
+    },
+
+    'DenseMatrix, Array | Matrix': function (m, b) {
+      return _denseForwardSubstitution(m, b)
+    },
+
+    'Array, Array | Matrix': function (a, b) {
+      const m = matrix(a)
+      const R = _denseForwardSubstitution(m, b)
+      return R.map(r => r.valueOf())
+    }
+  })
+
+  function _denseForwardSubstitution (m, b_) {
+    // the algorithm is derived from
+    // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93
+
+    // array of right-hand sides
+    const B = [solveValidation(m, b_, true)._data.map(e => e[0])]
+
+    const M = m._data
+    const rows = m._size[0]
+    const columns = m._size[1]
+
+    // loop columns
+    for (let i = 0; i < columns; i++) {
+      let L = B.length
+
+      // loop right-hand sides
+      for (let k = 0; k < L; k++) {
+        const b = B[k]
+
+        if (!equalScalar(M[i][i], 0)) {
+          // non-singular row
+
+          b[i] = divideScalar(b[i], M[i][i])
+
+          for (let j = i + 1; j < columns; j++) {
+            // b[j] -= b[i] * M[j,i]
+            b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i]))
+          }
+        } else if (!equalScalar(b[i], 0)) {
+          // singular row, nonzero RHS
+
+          if (k === 0) {
+            // There is no valid solution
+            return []
+          } else {
+            // This RHS is invalid but other solutions may still exist
+            B.splice(k, 1)
+            k -= 1
+            L -= 1
+          }
+        } else if (k === 0) {
+          // singular row, RHS is zero
+
+          const bNew = [...b]
+          bNew[i] = 1
+
+          for (let j = i + 1; j < columns; j++) {
+            bNew[j] = subtract(bNew[j], M[j][i])
+          }
+
+          B.push(bNew)
+        }
+      }
+    }
+
+    return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] }))
+  }
+
+  function _sparseForwardSubstitution (m, b_) {
+    // array of right-hand sides
+    const B = [solveValidation(m, b_, true)._data.map(e => e[0])]
+
+    const rows = m._size[0]
+    const columns = m._size[1]
+
+    const values = m._values
+    const index = m._index
+    const ptr = m._ptr
+
+    // loop columns
+    for (let i = 0; i < columns; i++) {
+      let L = B.length
+
+      // loop right-hand sides
+      for (let k = 0; k < L; k++) {
+        const b = B[k]
+
+        // values & indices (column i)
+        const iValues = []
+        const iIndices = []
+
+        // first & last indeces in column
+        const firstIndex = ptr[i]
+        const lastIndex = ptr[i + 1]
+
+        // find the value at [i, i]
+        let Mii = 0
+        for (let j = firstIndex; j < lastIndex; j++) {
+          const J = index[j]
+          // check row
+          if (J === i) {
+            Mii = values[j]
+          } else if (J > i) {
+            // store lower triangular
+            iValues.push(values[j])
+            iIndices.push(J)
+          }
+        }
+
+        if (!equalScalar(Mii, 0)) {
+          // non-singular row
+
+          b[i] = divideScalar(b[i], Mii)
+
+          for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) {
+            const J = iIndices[j]
+            b[J] = subtract(b[J], multiplyScalar(b[i], iValues[j]))
+          }
+        } else if (!equalScalar(b[i], 0)) {
+          // singular row, nonzero RHS
+
+          if (k === 0) {
+            // There is no valid solution
+            return []
+          } else {
+            // This RHS is invalid but other solutions may still exist
+            B.splice(k, 1)
+            k -= 1
+            L -= 1
+          }
+        } else if (k === 0) {
+          // singular row, RHS is zero
+
+          const bNew = [...b]
+          bNew[i] = 1
+
+          for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) {
+            const J = iIndices[j]
+            bNew[J] = subtract(bNew[J], iValues[j])
+          }
+
+          B.push(bNew)
+        }
+      }
+    }
+
+    return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] }))
+  }
+})
diff --git a/src/function/algebra/solver/lusolve.js b/src/function/algebra/solver/lusolve.js
index fe64bf00ae..8bebb86051 100644
--- a/src/function/algebra/solver/lusolve.js
+++ b/src/function/algebra/solver/lusolve.js
@@ -53,67 +53,57 @@ export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ type
   return typed(name, {
 
     'Array, Array | Matrix': function (a, b) {
-      // convert a to matrix
       a = matrix(a)
-      // matrix lup decomposition
       const d = lup(a)
-      // solve
       const x = _lusolve(d.L, d.U, d.p, null, b)
-      // convert result to array
       return x.valueOf()
     },
 
     'DenseMatrix, Array | Matrix': function (a, b) {
-      // matrix lup decomposition
       const d = lup(a)
-      // solve
       return _lusolve(d.L, d.U, d.p, null, b)
     },
 
     'SparseMatrix, Array | Matrix': function (a, b) {
-      // matrix lup decomposition
       const d = lup(a)
-      // solve
       return _lusolve(d.L, d.U, d.p, null, b)
     },
 
     'SparseMatrix, Array | Matrix, number, number': function (a, b, order, threshold) {
-      // matrix lu decomposition
       const d = slu(a, order, threshold)
-      // solve
       return _lusolve(d.L, d.U, d.p, d.q, b)
     },
 
     'Object, Array | Matrix': function (d, b) {
-      // solve
       return _lusolve(d.L, d.U, d.p, d.q, b)
     }
   })
 
   function _toMatrix (a) {
-    // check it is a matrix
     if (isMatrix(a)) { return a }
-    // check array
     if (isArray(a)) { return matrix(a) }
-    // throw
     throw new TypeError('Invalid Matrix LU decomposition')
   }
 
   function _lusolve (l, u, p, q, b) {
-    // verify L, U, P
+    // verify decomposition
     l = _toMatrix(l)
     u = _toMatrix(u)
-    // validate matrix and vector
-    b = solveValidation(l, b, false)
+
     // apply row permutations if needed (b is a DenseMatrix)
-    if (p) { b._data = csIpvec(p, b._data) }
+    if (p) {
+      b = solveValidation(l, b, true)
+      b._data = csIpvec(p, b._data)
+    }
+
     // use forward substitution to resolve L * y = b
     const y = lsolve(l, b)
     // use backward substitution to resolve U * x = y
     const x = usolve(u, y)
+
     // apply column permutations if needed (x is a DenseMatrix)
     if (q) { x._data = csIpvec(q, x._data) }
-    // return solution
+
     return x
   }
 })
diff --git a/src/function/algebra/solver/usolve.js b/src/function/algebra/solver/usolve.js
index 8767ba1753..87bc9c93b2 100644
--- a/src/function/algebra/solver/usolve.js
+++ b/src/function/algebra/solver/usolve.js
@@ -16,7 +16,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
   const solveValidation = createSolveValidation({ DenseMatrix })
 
   /**
-   * Solves the linear equation system by backward substitution. Matrix must be an upper triangular matrix.
+   * Finds one solution of a linear equation system by backward substitution. Matrix must be an upper triangular matrix. Throws an error if there's no solution.
    *
    * `U * x = b`
    *
@@ -32,7 +32,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *
    * See also:
    *
-   *    lup, slu, usolve, lusolve
+   *    usolveAll, lup, slu, usolve, lusolve
    *
    * @param {Matrix, Array} U       A N x N matrix or array (U)
    * @param {Matrix, Array} b       A column vector with the b values
@@ -42,67 +42,64 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
   return typed(name, {
 
     'SparseMatrix, Array | Matrix': function (m, b) {
-      // process matrix
       return _sparseBackwardSubstitution(m, b)
     },
 
     'DenseMatrix, Array | Matrix': function (m, b) {
-      // process matrix
       return _denseBackwardSubstitution(m, b)
     },
 
     'Array, Array | Matrix': function (a, b) {
-      // create dense matrix from array
       const m = matrix(a)
-      // use matrix implementation
       const r = _denseBackwardSubstitution(m, b)
-      // result
       return r.valueOf()
     }
   })
 
   function _denseBackwardSubstitution (m, b) {
-    // validate matrix and vector, return copy of column vector b
+    // make b into a column vector
     b = solveValidation(m, b, true)
-    // column vector data
+
     const bdata = b._data
-    // rows & columns
+
     const rows = m._size[0]
     const columns = m._size[1]
+
     // result
     const x = []
-    // arrays
-    const data = m._data
-    // backward solve m * x = b, loop columns (backwards)
+
+    const mdata = m._data
+    // loop columns backwards
     for (let j = columns - 1; j >= 0; j--) {
       // b[j]
       const bj = bdata[j][0] || 0
       // x[j]
       let xj
-      // backward substitution (outer product) avoids inner looping when bj === 0
+
       if (!equalScalar(bj, 0)) {
-        // value @ [j, j]
-        const vjj = data[j][j]
-        // check vjj
+        // value at [j, j]
+        const vjj = mdata[j][j]
+
         if (equalScalar(vjj, 0)) {
           // system cannot be solved
           throw new Error('Linear system cannot be solved since matrix is singular')
         }
-        // calculate xj
+
         xj = divideScalar(bj, vjj)
+
         // loop rows
         for (let i = j - 1; i >= 0; i--) {
           // update copy of b
-          bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, data[i][j]))]
+          bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))]
         }
       } else {
-        // zero value @ j
+        // zero value at j
         xj = 0
       }
       // update x
       x[j] = [xj]
     }
-    // return column vector
+
     return new DenseMatrix({
       data: x,
       size: [rows, 1]
@@ -110,71 +107,71 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed
   }
 
   function _sparseBackwardSubstitution (m, b) {
-    // validate matrix and vector, return copy of column vector b
+    // make b into a column vector
     b = solveValidation(m, b, true)
-    // column vector data
+
     const bdata = b._data
-    // rows & columns
+
     const rows = m._size[0]
     const columns = m._size[1]
-    // matrix arrays
+
     const values = m._values
     const index = m._index
     const ptr = m._ptr
-    // vars
-    let i, k
+
     // result
     const x = []
-    // backward solve m * x = b, loop columns (backwards)
+
+    // loop columns backwards
     for (let j = columns - 1; j >= 0; j--) {
-      // b[j]
       const bj = bdata[j][0] || 0
-      // backward substitution (outer product) avoids inner looping when bj === 0
+
       if (!equalScalar(bj, 0)) {
-        // value @ [j, j]
+        // non-degenerate row, find solution
+
         let vjj = 0
+
         // upper triangular matrix values & index (column j)
-        const jvalues = []
-        const jindex = []
+        const jValues = []
+        const jIndices = []
+
         // first & last indeces in column
-        const f = ptr[j]
-        let l = ptr[j + 1]
-        // values in column, find value @ [j, j], loop backwards
-        for (k = l - 1; k >= f; k--) {
-          // row
-          i = index[k]
-          // check row
+        const firstIndex = ptr[j]
+        const lastIndex = ptr[j + 1]
+
+        // values in column, find value at [j, j], loop backwards
+        for (let k = lastIndex - 1; k >= firstIndex; k--) {
+          const i = index[k]
+
+          // check row (rows are not sorted!)
           if (i === j) {
-            // update vjj
             vjj = values[k]
           } else if (i < j) {
             // store upper triangular
-            jvalues.push(values[k])
-            jindex.push(i)
+            jValues.push(values[k])
+            jIndices.push(i)
           }
         }
-        // at this point we must have a value @ [j, j]
+
+        // at this point we must have a value in vjj
         if (equalScalar(vjj, 0)) {
-          // system cannot be solved, there is no value @ [j, j]
           throw new Error('Linear system cannot be solved since matrix is singular')
         }
-        // calculate xj
+
         const xj = divideScalar(bj, vjj)
-        // loop upper triangular
-        for (k = 0, l = jindex.length; k < l; k++) {
-          // row
-          i = jindex[k]
-          // update copy of b
-          bdata[i] = [subtract(bdata[i][0], multiplyScalar(xj, jvalues[k]))]
+
+        for (let k = 0, lastIndex = jIndices.length; k < lastIndex; k++) {
+          const i = jIndices[k]
+          bdata[i] = [subtract(bdata[i][0], multiplyScalar(xj, jValues[k]))]
         }
-        // update x
+
         x[j] = [xj]
       } else {
-        // update x
+        // degenerate row, we can choose any value
         x[j] = [0]
       }
     }
-    // return vector
+
     return new DenseMatrix({
       data: x,
       size: [rows, 1]
diff --git a/src/function/algebra/solver/usolveAll.js b/src/function/algebra/solver/usolveAll.js
new file mode 100644
index 0000000000..9e1b9822fd
--- /dev/null
+++ b/src/function/algebra/solver/usolveAll.js
@@ -0,0 +1,199 @@
+import { factory } from '../../../utils/factory'
+import { createSolveValidation } from './utils/solveValidation'
+
+const name = 'usolveAll'
+const dependencies = [
+  'typed',
+  'matrix',
+  'divideScalar',
+  'multiplyScalar',
+  'subtract',
+  'equalScalar',
+  'DenseMatrix'
+]
+
+export const createUsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, divideScalar, multiplyScalar, subtract, equalScalar, DenseMatrix }) => {
+  const solveValidation = createSolveValidation({ DenseMatrix })
+
+  /**
+   * Finds all solutions of a linear equation system by backward substitution. Matrix must be an upper triangular matrix.
+   *
+   * `U * x = b`
+   *
+   * Syntax:
+   *
+   *    math.usolve(U, b)
+   *
+   * Examples:
+   *
+   *    const a = [[-2, 3], [2, 1]]
+   *    const b = [11, 9]
+   *    const x = usolve(a, b)  // [ [[8], [9]] ]
+   *
+   * See also:
+   *
+   *    usolve, lup, slu, usolve, lusolve
+   *
+   * @param {Matrix, Array} U       A N x N matrix or array (U)
+   * @param {Matrix, Array} b       A column vector with the b values
+   *
+   * @return {DenseMatrix[] | Array[]}  An array of affine-independent column vectors (x) that solve the linear system
+   */
+  return typed(name, {
+
+    'SparseMatrix, Array | Matrix': function (m, b) {
+      return _sparseBackwardSubstitution(m, b)
+    },
+
+    'DenseMatrix, Array | Matrix': function (m, b) {
+      return _denseBackwardSubstitution(m, b)
+    },
+
+    'Array, Array | Matrix': function (a, b) {
+      const m = matrix(a)
+      const R = _denseBackwardSubstitution(m, b)
+      return R.map(r => r.valueOf())
+    }
+  })
+
+  function _denseBackwardSubstitution (m, b_) {
+    // the algorithm is derived from
+    // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93
+
+    // array of right-hand sides
+    const B = [solveValidation(m, b_, true)._data.map(e => e[0])]
+
+    const M = m._data
+    const rows = m._size[0]
+    const columns = m._size[1]
+
+    // loop columns backwards
+    for (let i = columns - 1; i >= 0; i--) {
+      let L = B.length
+
+      // loop right-hand sides
+      for (let k = 0; k < L; k++) {
+        const b = B[k]
+
+        if (!equalScalar(M[i][i], 0)) {
+          // non-singular row
+
+          b[i] = divideScalar(b[i], M[i][i])
+
+          for (let j = i - 1; j >= 0; j--) {
+            // b[j] -= b[i] * M[j,i]
+            b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i]))
+          }
+        } else if (!equalScalar(b[i], 0)) {
+          // singular row, nonzero RHS
+
+          if (k === 0) {
+            // There is no valid solution
+            return []
+          } else {
+            // This RHS is invalid but other solutions may still exist
+            B.splice(k, 1)
+            k -= 1
+            L -= 1
+          }
+        } else if (k === 0) {
+          // singular row, RHS is zero
+
+          const bNew = [...b]
+          bNew[i] = 1
+
+          for (let j = i - 1; j >= 0; j--) {
+            bNew[j] = subtract(bNew[j], M[j][i])
+          }
+
+          B.push(bNew)
+        }
+      }
+    }
+
+    return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] }))
+  }
+
+  function _sparseBackwardSubstitution (m, b_) {
+    // array of right-hand sides
+    const B = [solveValidation(m, b_, true)._data.map(e => e[0])]
+
+    const rows = m._size[0]
+    const columns = m._size[1]
+
+    const values = m._values
+    const index = m._index
+    const ptr = m._ptr
+
+    // loop columns backwards
+    for (let i = columns - 1; i >= 0; i--) {
+      let L = B.length
+
+      // loop right-hand sides
+      for (let k = 0; k < L; k++) {
+        const b = B[k]
+
+        // values & indices (column i)
+        const iValues = []
+        const iIndices = []
+
+        // first & last indeces in column
+        const firstIndex = ptr[i]
+        const lastIndex = ptr[i + 1]
+
+        // find the value at [i, i]
+        let Mii = 0
+        for (let j = lastIndex - 1; j >= firstIndex; j--) {
+          const J = index[j]
+          // check row
+          if (J === i) {
+            Mii = values[j]
+          } else if (J < i) {
+            // store upper triangular
+            iValues.push(values[j])
+            iIndices.push(J)
+          }
+        }
+
+        if (!equalScalar(Mii, 0)) {
+          // non-singular row
+
+          b[i] = divideScalar(b[i], Mii)
+
+          // loop upper triangular
+          for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) {
+            const J = iIndices[j]
+            b[J] = subtract(b[J], multiplyScalar(b[i], iValues[j]))
+          }
+        } else if (!equalScalar(b[i], 0)) {
+          // singular row, nonzero RHS
+
+          if (k === 0) {
+            // There is no valid solution
+            return []
+          } else {
+            // This RHS is invalid but other solutions may still exist
+            B.splice(k, 1)
+            k -= 1
+            L -= 1
+          }
+        } else if (k === 0) {
+          // singular row, RHS is zero
+
+          const bNew = [...b]
+          bNew[i] = 1
+
+          // loop upper triangular
+          for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) {
+            const J = iIndices[j]
+            bNew[J] = subtract(bNew[J], iValues[j])
+          }
+
+          B.push(bNew)
+        }
+      }
+    }
+
+    return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] }))
+  }
+})
diff --git a/src/function/algebra/solver/utils/solveValidation.js b/src/function/algebra/solver/utils/solveValidation.js
index 056a8f9db0..eb6dd66965 100644
--- a/src/function/algebra/solver/utils/solveValidation.js
+++ b/src/function/algebra/solver/utils/solveValidation.js
@@ -1,4 +1,4 @@
-import { isArray, isDenseMatrix, isMatrix } from '../../../../utils/is'
+import { isArray, isMatrix, isDenseMatrix, isSparseMatrix } from '../../../../utils/is'
 import { arraySize } from '../../../../utils/array'
 import { format } from '../../../../utils/string'
 
@@ -13,131 +13,123 @@ export function createSolveValidation ({ DenseMatrix }) {
    * @return {DenseMatrix}        Dense column vector b
    */
   return function solveValidation (m, b, copy) {
-    // matrix size
-    const size = m.size()
-    // validate matrix dimensions
-    if (size.length !== 2) { throw new RangeError('Matrix must be two dimensional (size: ' + format(size) + ')') }
-    // rows & columns
-    const rows = size[0]
-    const columns = size[1]
-    // validate rows & columns
-    if (rows !== columns) { throw new RangeError('Matrix must be square (size: ' + format(size) + ')') }
-    // vars
-    let data, i, bdata
-    // check b is matrix
+    const mSize = m.size()
+
+    if (mSize.length !== 2) {
+      throw new RangeError('Matrix must be two dimensional (size: ' + format(mSize) + ')')
+    }
+
+    const rows = mSize[0]
+    const columns = mSize[1]
+
+    if (rows !== columns) {
+      throw new RangeError('Matrix must be square (size: ' + format(mSize) + ')')
+    }
+
+    let data = []
+
     if (isMatrix(b)) {
-      // matrix size
-      const msize = b.size()
-      // vector
-      if (msize.length === 1) {
-        // check vector length
-        if (msize[0] !== rows) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') }
-        // create data array
-        data = []
-        // matrix data (DenseMatrix)
-        bdata = b._data
-        // loop b data
-        for (i = 0; i < rows; i++) {
-          // row array
+      const bSize = b.size()
+      const bdata = b._data
+
+      // 1-dim vector
+      if (bSize.length === 1) {
+        if (bSize[0] !== rows) {
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.')
+        }
+
+        for (let i = 0; i < rows; i++) {
           data[i] = [bdata[i]]
         }
-        // return Dense Matrix
+
         return new DenseMatrix({
           data: data,
           size: [rows, 1],
           datatype: b._datatype
         })
       }
-      // two dimensions
-      if (msize.length === 2) {
-        // array must be a column vector
-        if (msize[0] !== rows || msize[1] !== 1) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') }
-        // check matrix type
+
+      // 2-dim column
+      if (bSize.length === 2) {
+        if (bSize[0] !== rows || bSize[1] !== 1) {
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.')
+        }
+
         if (isDenseMatrix(b)) {
-          // check a copy is needed
           if (copy) {
-            // create data array
             data = []
-            // matrix data (DenseMatrix)
-            bdata = b._data
-            // loop b data
-            for (i = 0; i < rows; i++) {
-              // row array
+
+            for (let i = 0; i < rows; i++) {
               data[i] = [bdata[i][0]]
             }
-            // return Dense Matrix
+
             return new DenseMatrix({
               data: data,
               size: [rows, 1],
               datatype: b._datatype
             })
           }
-          // b is already a column vector
+
           return b
         }
-        // create data array
-        data = []
-        for (i = 0; i < rows; i++) { data[i] = [0] }
-        // sparse matrix arrays
-        const values = b._values
-        const index = b._index
-        const ptr = b._ptr
-        // loop values in column 0
-        for (let k1 = ptr[1], k = ptr[0]; k < k1; k++) {
-          // row
-          i = index[k]
-          // add to data
-          data[i][0] = values[k]
+
+        if (isSparseMatrix(b)) {
+          for (let i = 0; i < rows; i++) { data[i] = [0] }
+
+          const values = b._values
+          const index = b._index
+          const ptr = b._ptr
+
+          for (let k1 = ptr[1], k = ptr[0]; k < k1; k++) {
+            const i = index[k]
+            data[i][0] = values[k]
+          }
+
+          return new DenseMatrix({
+            data: data,
+            size: [rows, 1],
+            datatype: b._datatype
+          })
         }
-        // return Dense Matrix
-        return new DenseMatrix({
-          data: data,
-          size: [rows, 1],
-          datatype: b._datatype
-        })
       }
-      // throw error
-      throw new RangeError('Dimension mismatch. Matrix columns must match vector length.')
+
+      throw new RangeError('Dimension mismatch. The right side has to be either 1- or 2-dimensional vector.')
     }
-    // check b is array
+
     if (isArray(b)) {
-      // size
-      const asize = arraySize(b)
-      // check matrix dimensions, vector
-      if (asize.length === 1) {
-        // check vector length
-        if (asize[0] !== rows) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') }
-        // create data array
-        data = []
-        // loop b
-        for (i = 0; i < rows; i++) {
-          // row array
+      const bsize = arraySize(b)
+
+      if (bsize.length === 1) {
+        if (bsize[0] !== rows) {
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.')
+        }
+
+        for (let i = 0; i < rows; i++) {
           data[i] = [b[i]]
         }
-        // return Dense Matrix
+
         return new DenseMatrix({
           data: data,
           size: [rows, 1]
         })
       }
-      if (asize.length === 2) {
-        // array must be a column vector
-        if (asize[0] !== rows || asize[1] !== 1) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') }
-        // create data array
-        data = []
-        // loop b data
-        for (i = 0; i < rows; i++) {
-          // row array
+
+      if (bsize.length === 2) {
+        if (bsize[0] !== rows || bsize[1] !== 1) {
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.')
+        }
+
+        for (let i = 0; i < rows; i++) {
           data[i] = [b[i][0]]
         }
-        // return Dense Matrix
+
         return new DenseMatrix({
           data: data,
           size: [rows, 1]
         })
       }
-      // throw error
-      throw new RangeError('Dimension mismatch. Matrix columns must match vector length.')
+
+      throw new RangeError('Dimension mismatch. The right side has to be either 1- or 2-dimensional vector.')
     }
   }
 }
diff --git a/test/unit-tests/function/algebra/solver/lsolveAll.test.js b/test/unit-tests/function/algebra/solver/lsolveAll.test.js
new file mode 100644
index 0000000000..6df922ba84
--- /dev/null
+++ b/test/unit-tests/function/algebra/solver/lsolveAll.test.js
@@ -0,0 +1,155 @@
+// test lsolveAll
+import assert from 'assert'
+
+import approx from '../../../../../tools/approx'
+import math from '../../../../../src/bundleAny'
+
+describe('lsolveAll', function () {
+  it('should solve linear system 4 x 4, arrays', function () {
+    const m =
+        [
+          [1, 0, 0, 0],
+          [1, 1, 0, 0],
+          [1, 1, 1, 0],
+          [1, 1, 1, 1]
+        ]
+    const b = [1, 2, 3, 4]
+
+    const x = math.lsolveAll(m, b)
+
+    approx.deepEqual(x, [[[1], [1], [1], [1]]])
+  })
+
+  it('should solve linear system 4 x 4, array and column array', function () {
+    const m =
+        [
+          [1, 0, 0, 0],
+          [1, 1, 0, 0],
+          [1, 1, 1, 0],
+          [1, 1, 1, 1]
+        ]
+    const b = [
+      [1],
+      [2],
+      [3],
+      [4]
+    ]
+    const x = math.lsolveAll(m, b)
+
+    approx.deepEqual(x, [[[1], [1], [1], [1]]])
+  })
+
+  it('should solve linear system 4 x 4, matrices', function () {
+    const m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [1, 1, 0, 0],
+        [1, 1, 1, 0],
+        [1, 1, 1, 1]
+      ])
+    const b = math.matrix([1, 2, 3, 4])
+
+    const x = math.lsolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])])
+  })
+
+  it('should solve linear system 4 x 4, sparse matrices', function () {
+    const m = math.sparse(
+      [
+        [1, 0, 0, 0],
+        [1, 1, 0, 0],
+        [1, 1, 1, 0],
+        [1, 1, 1, 1]
+      ])
+    const b = math.matrix([[1], [2], [3], [4]], 'sparse')
+
+    const x = math.lsolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])])
+  })
+
+  it('should solve linear system 4 x 4, matrix and column matrix', function () {
+    const m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [1, 1, 0, 0],
+        [1, 1, 1, 0],
+        [1, 1, 1, 1]
+      ])
+    const b = math.matrix([
+      [1],
+      [2],
+      [3],
+      [4]
+    ])
+
+    const x = math.lsolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])])
+  })
+
+  it('should solve linear system 4 x 4, sparse matrix and column matrix', function () {
+    const m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [1, 1, 0, 0],
+        [1, 1, 1, 0],
+        [1, 1, 1, 1]
+      ], 'sparse')
+    const b = math.matrix([
+      [1],
+      [2],
+      [3],
+      [4]
+    ], 'sparse')
+
+    const x = math.lsolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])])
+  })
+
+  it('should return an empty array when there is no solution', function () {
+    assert.deepStrictEqual([], math.lsolveAll([[1, 1], [0, 0]], [1, 1]))
+    assert.deepStrictEqual([], math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]))
+    assert.deepStrictEqual([], math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]))
+  })
+
+  it('should solve systems with singular dense matrices', function () {
+    approx.deepEqual(
+      math.lsolveAll([[2, 0, 0], [1, 0, 0], [-1, 1, 1]], [4, 2, 1]),
+      [[[2], [0], [3]], [[2], [1], [2]]]
+    )
+
+    approx.deepEqual(
+      math.lsolveAll([[0, 0, 0], [1, 1, 0], [2, 1, 0]], [0, 2, 2]),
+      [[[0], [2], [0]], [[0], [2], [1]]]
+    )
+
+    approx.deepEqual(
+      math.lsolveAll([[0, 0, 0], [1, 1, 0], [1, 1, 0]], [0, 2, 2]),
+      [[[0], [2], [0]], [[1], [1], [0]], [[0], [2], [1]]]
+    )
+  })
+
+  it('should solve systems with singular sparse matrices', function () {
+    approx.deepEqual(
+      math.lsolveAll(math.matrix([[2, 0, 0], [1, 0, 0], [-1, 1, 1]], 'sparse'), [4, 2, 1]),
+      [math.matrix([[2], [0], [3]]), math.matrix([[2], [1], [2]])]
+    )
+
+    approx.deepEqual(
+      math.lsolveAll(math.matrix([[0, 0, 0], [1, 1, 0], [2, 1, 0]], 'sparse'), [0, 2, 2]),
+      [math.matrix([[0], [2], [0]]), math.matrix([[0], [2], [1]])]
+    )
+
+    approx.deepEqual(
+      math.lsolveAll(math.matrix([[0, 0, 0], [1, 1, 0], [1, 1, 0]], 'sparse'), [0, 2, 2]),
+      [math.matrix([[0], [2], [0]]), math.matrix([[1], [1], [0]]), math.matrix([[0], [2], [1]])]
+    )
+  })
+})
diff --git a/test/unit-tests/function/algebra/solver/usolveAll.test.js b/test/unit-tests/function/algebra/solver/usolveAll.test.js
new file mode 100644
index 0000000000..f7e3453146
--- /dev/null
+++ b/test/unit-tests/function/algebra/solver/usolveAll.test.js
@@ -0,0 +1,155 @@
+// test usolveAll
+import assert from 'assert'
+
+import approx from '../../../../../tools/approx'
+import math from '../../../../../src/bundleAny'
+
+describe('usolveAll', function () {
+  it('should solve linear system 4 x 4, arrays', function () {
+    const m =
+        [
+          [1, 1, 1, 1],
+          [0, 1, 1, 1],
+          [0, 0, 1, 1],
+          [0, 0, 0, 1]
+        ]
+    const b = [1, 2, 3, 4]
+
+    const x = math.usolveAll(m, b)
+
+    approx.deepEqual(x, [[[-1], [-1], [-1], [4]]])
+  })
+
+  it('should solve linear system 4 x 4, array and column array', function () {
+    const m =
+        [
+          [1, 1, 1, 1],
+          [0, 1, 1, 1],
+          [0, 0, 1, 1],
+          [0, 0, 0, 1]
+        ]
+    const b = [
+      [1],
+      [2],
+      [3],
+      [4]
+    ]
+    const x = math.usolveAll(m, b)
+
+    approx.deepEqual(x, [[[-1], [-1], [-1], [4]]])
+  })
+
+  it('should solve linear system 4 x 4, matrices', function () {
+    const m = math.matrix(
+      [
+        [1, 1, 1, 1],
+        [0, 1, 1, 1],
+        [0, 0, 1, 1],
+        [0, 0, 0, 1]
+      ])
+    const b = math.matrix([1, 2, 3, 4])
+
+    const x = math.usolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])])
+  })
+
+  it('should solve linear system 4 x 4, sparse matrices', function () {
+    const m = math.sparse(
+      [
+        [1, 1, 1, 1],
+        [0, 1, 1, 1],
+        [0, 0, 1, 1],
+        [0, 0, 0, 1]
+      ])
+    const b = math.matrix([[1], [2], [3], [4]], 'sparse')
+
+    const x = math.usolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])])
+  })
+
+  it('should solve linear system 4 x 4, matrix and column matrix', function () {
+    const m = math.matrix(
+      [
+        [1, 1, 1, 1],
+        [0, 1, 1, 1],
+        [0, 0, 1, 1],
+        [0, 0, 0, 1]
+      ])
+    const b = math.matrix([
+      [1],
+      [2],
+      [3],
+      [4]
+    ])
+
+    const x = math.usolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])])
+  })
+
+  it('should solve linear system 4 x 4, sparse matrix and column matrix', function () {
+    const m = math.matrix(
+      [
+        [1, 1, 1, 1],
+        [0, 1, 1, 1],
+        [0, 0, 1, 1],
+        [0, 0, 0, 1]
+      ], 'sparse')
+    const b = math.matrix([
+      [1],
+      [2],
+      [3],
+      [4]
+    ], 'sparse')
+
+    const x = math.usolveAll(m, b)
+
+    assert(x[0] instanceof math.Matrix)
+    approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])])
+  })
+
+  it('should return an empty array when there is no solution', function () {
+    assert.deepStrictEqual([], math.usolveAll([[1, 1], [0, 0]], [1, 1]))
+    assert.deepStrictEqual([], math.usolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]))
+    assert.deepStrictEqual([], math.usolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]))
+  })
+
+  it('should solve systems with singular dense matrices', function () {
+    approx.deepEqual(
+      math.usolveAll([[1, 1, -1], [0, 0, 1], [0, 0, 2]], [1, 2, 4]),
+      [[[3], [0], [2]], [[2], [1], [2]]]
+    )
+
+    approx.deepEqual(
+      math.usolveAll([[0, 1, 2], [0, 1, 1], [0, 0, 0]], [2, 2, 0]),
+      [[[0], [2], [0]], [[1], [2], [0]]]
+    )
+
+    approx.deepEqual(
+      math.usolveAll([[0, 1, 1], [0, 1, 1], [0, 0, 0]], [2, 2, 0]),
+      [[[0], [2], [0]], [[0], [1], [1]], [[1], [2], [0]]]
+    )
+  })
+
+  it('should solve systems with singular sparse matrices', function () {
+    approx.deepEqual(
+      math.usolveAll(math.matrix([[1, 1, -1], [0, 0, 1], [0, 0, 2]], 'sparse'), [1, 2, 4]),
+      [math.matrix([[3], [0], [2]]), math.matrix([[2], [1], [2]])]
+    )
+
+    approx.deepEqual(
+      math.usolveAll(math.matrix([[0, 1, 2], [0, 1, 1], [0, 0, 0]], 'sparse'), [2, 2, 0]),
+      [math.matrix([[0], [2], [0]]), math.matrix([[1], [2], [0]])]
+    )
+
+    approx.deepEqual(
+      math.usolveAll(math.matrix([[0, 1, 1], [0, 1, 1], [0, 0, 0]], 'sparse'), [2, 2, 0]),
+      [math.matrix([[0], [2], [0]]), math.matrix([[0], [1], [1]]), math.matrix([[1], [2], [0]])]
+    )
+  })
+})