diff --git a/AUTHORS b/AUTHORS
index 88eba60aff..7ca6b9cece 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -263,6 +263,7 @@ witer33 <34513257+witer33@users.noreply.github.com>
 Hudsxn <jhudson98new@gmail.com>
 Christian Stussak <stussak@mfo.de>
 aitee <1228639+aitee@users.noreply.github.com>
+Kip Robinson <91914404+kiprobinsonknack@users.noreply.github.com>
 Delaney Sylvans <delaneysylvans@gmail.com>
 Rani D. <73716554+ranidam@users.noreply.github.com>
 Don McCurdy <dm@donmccurdy.com>
diff --git a/HISTORY.md b/HISTORY.md
index 43afcc5430..8582126b3b 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -1,5 +1,22 @@
 # History
 
+# Unpublished changes scheduled for v15
+
+!!! BE CAREFUL: BREAKING CHANGES !!!
+
+- Fix: #1753 Correct dimensionality of Kronecker product on vectors (and
+  extend to arbitrary dimension) (#3455). Thanks @Delaney.
+- Fix: #3501 parse `%` as unary only when not followed by a term (#3505).
+  Thanks @gwhitney.
+- Fix #3421 require a space or delimiter after hex, bin, and oct values (#3463).
+- Feat: #3349 Decouple precedence of unary percentage operator and binary
+  modulus operator (that both use symbol `%`), and raise the former (#3432).
+  Thanks @kiprobinsonknack.
+- Feat: #1753 enhance Kronecker product to handle arbitrary dimension (#3461).
+  Thanks @gwhitney.
+- Feat: #2344 matrix subset according to the type of input (#3485).
+  Thanks @dvd101x.
+
 # 2025-09-26, 14.8.1
 
 - Fix: #3538 `config` printing a warning when using `{ number: 'bigint' }` 
diff --git a/docs/core/configuration.md b/docs/core/configuration.md
index b72fdcc81e..95f93329c1 100644
--- a/docs/core/configuration.md
+++ b/docs/core/configuration.md
@@ -15,7 +15,8 @@ const config = {
   numberFallback: 'number',
   precision: 64,
   predictable: false,
-  randomSeed: null
+  randomSeed: null,
+  legacySubset: false
 }
 const math = create(all, config)
 
@@ -84,6 +85,8 @@ The following configuration options are available:
 
 - `randomSeed`. Set this option to seed pseudo random number generation, making it deterministic. The pseudo random number generator is reset with the seed provided each time this option is set. For example, setting it to `'a'` will cause `math.random()` to return `0.43449421599986604` upon the first call after setting the option every time. Set to `null` to seed the pseudo random number generator with a random seed. Default value is `null`.
 
+- `legacySubset`. When set to `true`, the `subset` function behaves as in earlier versions of math.js: retrieving a subset where the index size contains only one element, returns the value itself. When set to `false` (default), `subset` eliminates dimensions from the result only if the index dimension is a scalar; if the index dimension is a range, matrix, or array, the dimensions are preserved. This option is helpful for maintaining compatibility with legacy code that depends on the previous behavior.
+
 
 ## Examples
 
diff --git a/docs/datatypes/matrices.md b/docs/datatypes/matrices.md
index 4f888d6b4b..b362a4ead6 100644
--- a/docs/datatypes/matrices.md
+++ b/docs/datatypes/matrices.md
@@ -270,13 +270,34 @@ in the matrix, and if not, a subset of the matrix will be returned.
 
 A subset can be defined using an `Index`. An `Index` contains a single value
 or a set of values for each dimension of a matrix. An `Index` can be
-created using the function `index`. When getting a single value from a matrix,
-`subset` will return the value itself instead of a matrix containing just this 
-value.
+created using the function `index`. The way `subset` returns results depends on how you specify indices for each dimension:
 
-The function `subset` normally returns a subset, but when getting or setting a
-single value in a matrix, the value itself is returned.
+- If you use a scalar (single number) as an index for a dimension, that dimension is removed from the result.
+- If you use an array, matrix or range (even with just one element) as an index, that dimension is preserved in the result.
 
+This means that scalar indices eliminate dimensions, while array, matrix or range indices retain them. See the section [Migrate to v15](#migrate-to-v15) for more details and examples of this behavior.
+
+For example:
+
+```js
+const m = [
+  [10, 11, 12],
+  [20, 21, 22]
+]
+
+// Scalar index eliminates the dimension:
+math.subset(m, math.index(1, 2))           // 22 (both dimensions indexed by scalars, result is a value)
+math.subset(m, math.index(1, [2]))         // [22] (row dimension eliminated, column dimension preserved as array)
+math.subset(m, math.index([1], 2))         // [22] (column dimension eliminated, row dimension preserved as array)
+math.subset(m, math.index([1], [2]))       // [[22]] (both dimensions preserved as arrays)
+
+math.config({legacySubset: true}) // switch to legacy behavior
+math.subset(m, math.index(1, 2))           // 22
+math.subset(m, math.index(1, [2]))         // 22
+math.subset(m, math.index([1], 2))         // 22
+math.subset(m, math.index([1], [2]))       // 22
+
+```
 
 Matrix indexes in math.js are zero-based, like most programming languages
 including JavaScript itself. Note that mathematical applications like Matlab 
@@ -296,7 +317,7 @@ math.subset(a, math.index([2, 3]))            // Array, [2, 3]
 math.subset(a, math.index(math.range(0,4)))   // Array, [0, 1, 2, 3]
 math.subset(b, math.index(1, 0))              // 2
 math.subset(b, math.index(1, [0, 1]))         // Array, [2, 3]
-math.subset(b, math.index([0, 1], 0))         // Matrix, [[0], [2]]
+math.subset(b, math.index([0, 1], [0]))       // Matrix, [[0], [2]]
 
 // get a subset
 d.subset(math.index([1, 2], [0, 1]))          // Matrix, [[3, 4], [6, 7]]
@@ -313,6 +334,34 @@ e.resize([2, 3], 0)                           // Matrix, [[0, 0, 0], [0, 0, 0]]
 e.subset(math.index(1, 2), 5)                 // Matrix, [[0, 0, 0], [0, 0, 5]]
 ```
 
+## Migrate indexing behavior to mathjs v15
+
+With the release of math.js v15, the behavior of `subset` when indexing matrices and arrays has changed. If your code relies on the previous behavior (where indexing with an array or matrix of size 1 would always return the value itself), you may need to update your code or enable legacy mode.
+
+To maintain the old indexing behavior without need for any code changes, use the configuration option `legacySubset`:
+
+```js
+math.config({ legacySubset: true })
+```
+
+To migrate your code, you'll have to change all matrix indexes from the old index notation to the new index notation. Basically: scalar indexes have to be wrapped in array brackets if you want an array as output. Here some examples:
+
+```js
+const m = math.matrix([[1, 2, 3], [4, 5, 6]])
+```
+
+| v14 code                                     | v15 equivalent code                       | Result             |
+|----------------------------------------------|-------------------------------------------|--------------------|
+| `math.subset(m, math.index([0, 1], [1, 2]))` | No change needed                          | `[[2, 3], [5, 6]]` |
+| `math.subset(m, math.index(1, [1, 2]))`      | `math.subset(m, math.index([1], [1, 2]))` | `[[5, 6]]`         |
+| `math.subset(m, math.index([0, 1], 2))`      | `math.subset(m, math.index([0, 1], [2]))` | `[[3], [6]]`       |
+| `math.subset(m, math.index(1, 2))`           | No change needed                          | 6                  |
+
+
+> **Tip:**  
+> If you want to get a scalar value, use scalar indices.  
+> If you want to preserve dimensions, use array, matrix or range indices.
+
 ## Getting and setting a value in a matrix
 
 There are two methods available on matrices that allow to get or set a single 
diff --git a/docs/expressions/syntax.md b/docs/expressions/syntax.md
index ad8bd3c92d..f435773bd4 100644
--- a/docs/expressions/syntax.md
+++ b/docs/expressions/syntax.md
@@ -118,8 +118,9 @@ Operators                         | Description
 `??`                              | Nullish coalescing
 `^`, `.^`                         | Exponentiation
 `+`, `-`, `~`, `not`              | Unary plus, unary minus, bitwise not, logical not
+`%`                               | Unary percentage
 See section below                 | Implicit multiplication
-`*`, `/`, `.*`, `./`,`%`, `mod`   | Multiply, divide , percentage, modulus
+`*`, `/`, `.*`, `./`,`%`, `mod`   | Multiply, divide, modulus
 `+`, `-`                          | Add, subtract
 `:`                               | Range
 `to`, `in`                        | Unit conversion
@@ -138,8 +139,8 @@ See section below                 | Implicit multiplication
 `\n`, `;`                         | Statement separators
 
 Lazy evaluation is used where logically possible for bitwise and logical
-operators. In the following example, the value of `x` will not even be 
-evaluated because it cannot effect the final result: 
+operators. In the following example, the sub-expression `x` will not even be
+evaluated because it cannot affect the final result:
 ```js
 math.evaluate('false and x')        // false, no matter what x equals
 ```
@@ -677,13 +678,15 @@ at 1, when the end is undefined, the range will end at the end of the matrix.
 
 There is a context variable `end` available as well to denote the end of the
 matrix. This variable cannot be used in multiple nested indices. In that case,
-`end` will be resolved as the end of the innermost matrix. To solve this, 
+`end` will be resolved as the end of the innermost matrix. To solve this,
 resolving of the nested index needs to be split in two separate operations.
 
 *IMPORTANT: matrix indexes and ranges work differently from the math.js indexes
 in JavaScript: They are one-based with an included upper-bound, similar to most
 math applications.*
 
+*IMPORTANT: The matrix indexing behavior has been changed in mathjs `v15`, see section [Migrate matrix indexing to mathjs v15](#migrate-matrix-indexing-to-mathjs-v15).*
+
 ```js
 parser = math.parser()
 
@@ -701,10 +704,39 @@ parser.evaluate('d = a * b')              // Matrix, [[19, 22], [43, 50]]
 
 // retrieve a subset of a matrix
 parser.evaluate('d[2, 1]')                // 43
-parser.evaluate('d[2, 1:end]')            // Matrix, [[43, 50]]
+parser.evaluate('d[2, 1:end]')            // Matrix, [43, 50]
 parser.evaluate('c[end - 1 : -1 : 2]')    // Matrix, [8, 7, 6]
 ```
 
+#### Migrate matrix indexing to mathjs v15
+
+With mathjs v15, matrix indexing has changed to be more consistent and predictable. In v14, using a scalar index would sometimes reduce the dimensionality of the result. In v15, if you want to preserve dimensions, use array, matrix, or range indices. If you want a scalar value, use scalar indices.
+
+To maintain the old indexing behavior without need for any code changes, use the configuration option `legacySubset`:
+
+```js
+math.config({ legacySubset: true })
+```
+
+To migrate your code, you'll have to change all matrix indexes from the old index notation to the new index notation. Basically: scalar indexes have to be wrapped in array brackets if you want an array as output. Here some examples:
+
+```js
+parser = math.parser()
+parser.evaluate('m = [1, 2, 3; 4, 5, 6]')
+```
+
+| v14 code                | v15 equivalent code           | Result             |
+|-------------------------|-------------------------------|--------------------|
+| `m[1:2, 2:3]`           | No change needed              | `[[2, 3], [5, 6]]` |
+| `m[2, 2:3]`             | `m[[2], 2:3]`                 | `[[5, 6]]`         |
+| `m[1:2, 3]`             | `m[1:2, [3]]`                 | `[[3], [6]]`       |
+| `m[2, 3]`               | No change needed              | `6`                |
+
+> **Tip:**  
+> If you want to always get a scalar value, use scalar indices.  
+> If you want to preserve dimensions, use array, matrix or range indices.
+
+
 ## Objects
 
 Objects in math.js work the same as in languages like JavaScript and Python.
diff --git a/src/core/config.js b/src/core/config.js
index 005bba30f1..f777ae6181 100644
--- a/src/core/config.js
+++ b/src/core/config.js
@@ -28,5 +28,10 @@ export const DEFAULT_CONFIG = {
 
   // random seed for seeded pseudo random number generation
   // null = randomly seed
-  randomSeed: null
+  randomSeed: null,
+
+  // legacy behavior for matrix subset. When true, the subset function
+  // returns a matrix or array with the same size as the index (except for scalars).
+  // When false, it returns a matrix or array with a size depending on the type of index.
+  legacySubset: false
 }
diff --git a/src/core/function/config.js b/src/core/function/config.js
index e246de46f1..5be65558f6 100644
--- a/src/core/function/config.js
+++ b/src/core/function/config.js
@@ -60,6 +60,11 @@ export function configFactory (config, emit) {
         delete optionsFix.epsilon
         return _config(optionsFix)
       }
+
+      if (options.legacySubset === true) {
+        // this if is only for backwards compatibility, it can be removed in the future.
+        console.warn('Warning: The configuration option "legacySubset" is for compatibility only and might be deprecated in the future.')
+      }
       const prev = clone(config)
 
       // validate some of the options
diff --git a/src/expression/parse.js b/src/expression/parse.js
index 1b1ce48929..d3564cc76e 100644
--- a/src/expression/parse.js
+++ b/src/expression/parse.js
@@ -348,7 +348,10 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
         next(state)
         state.token += currentCharacter(state)
         next(state)
-        while (parse.isHexDigit(currentCharacter(state))) {
+        while (
+          parse.isAlpha(currentCharacter(state), prevCharacter(state), nextCharacter(state)) ||
+          parse.isDigit(currentCharacter(state))
+        ) {
           state.token += currentCharacter(state)
           next(state)
         }
@@ -357,7 +360,10 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
           state.token += '.'
           next(state)
           // get the digits after the radix
-          while (parse.isHexDigit(currentCharacter(state))) {
+          while (
+            parse.isAlpha(currentCharacter(state), prevCharacter(state), nextCharacter(state)) ||
+            parse.isDigit(currentCharacter(state))
+          ) {
             state.token += currentCharacter(state)
             next(state)
           }
@@ -576,17 +582,6 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
     return (c >= '0' && c <= '9')
   }
 
-  /**
-   * checks if the given char c is a hex digit
-   * @param {string} c   a string with one character
-   * @return {boolean}
-   */
-  parse.isHexDigit = function isHexDigit (c) {
-    return ((c >= '0' && c <= '9') ||
-            (c >= 'a' && c <= 'f') ||
-            (c >= 'A' && c <= 'F'))
-  }
-
   /**
    * Start of the parse levels below, in order of precedence
    * @return {Node} node
@@ -1003,7 +998,7 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
   function parseAddSubtract (state) {
     let node, name, fn, params
 
-    node = parseMultiplyDivideModulusPercentage(state)
+    node = parseMultiplyDivideModulus(state)
 
     const operators = {
       '+': 'add',
@@ -1014,7 +1009,7 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
       fn = operators[name]
 
       getTokenSkipNewline(state)
-      const rightNode = parseMultiplyDivideModulusPercentage(state)
+      const rightNode = parseMultiplyDivideModulus(state)
       if (rightNode.isPercentage) {
         params = [node, new OperatorNode('*', 'multiply', [node, rightNode])]
       } else {
@@ -1027,11 +1022,11 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
   }
 
   /**
-   * multiply, divide, modulus, percentage
+   * multiply, divide, modulus
    * @return {Node} node
    * @private
    */
-  function parseMultiplyDivideModulusPercentage (state) {
+  function parseMultiplyDivideModulus (state) {
     let node, last, name, fn
 
     node = parseImplicitMultiplication(state)
@@ -1051,25 +1046,9 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
         // explicit operators
         name = state.token
         fn = operators[name]
-
         getTokenSkipNewline(state)
-
-        if (name === '%' && state.tokenType === TOKENTYPE.DELIMITER && state.token !== '(') {
-          // If the expression contains only %, then treat that as /100
-          if (state.token !== '' && operators[state.token]) {
-            const left = new OperatorNode('/', 'divide', [node, new ConstantNode(100)], false, true)
-            name = state.token
-            fn = operators[name]
-            getTokenSkipNewline(state)
-            last = parseImplicitMultiplication(state)
-
-            node = new OperatorNode(name, fn, [left, last])
-          } else { node = new OperatorNode('/', 'divide', [node, new ConstantNode(100)], false, true) }
-          // return node
-        } else {
-          last = parseImplicitMultiplication(state)
-          node = new OperatorNode(name, fn, [node, last])
-        }
+        last = parseImplicitMultiplication(state)
+        node = new OperatorNode(name, fn, [node, last])
       } else {
         break
       }
@@ -1122,7 +1101,7 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
    * @private
    */
   function parseRule2 (state) {
-    let node = parseUnary(state)
+    let node = parseUnaryPercentage(state)
     let last = node
     const tokenStates = []
 
@@ -1145,7 +1124,7 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
             // Rewind once and build the "number / number" node; the symbol will be consumed later
             Object.assign(state, tokenStates.pop())
             tokenStates.pop()
-            last = parseUnary(state)
+            last = parseUnaryPercentage(state)
             node = new OperatorNode('/', 'divide', [node, last])
           } else {
             // Not a match, so rewind
@@ -1166,6 +1145,36 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
     return node
   }
 
+  /**
+   * Unary percentage operator (treated as `value / 100`)
+   * @return {Node} node
+   * @private
+   */
+  function parseUnaryPercentage (state) {
+    let node = parseUnary(state)
+
+    if (state.token === '%') {
+      const previousState = Object.assign({}, state)
+      getTokenSkipNewline(state)
+      // We need to decide if this is a unary percentage % or binary modulo %
+      // So we attempt to parse a unary expression at this point.
+      // If it fails, then the only possibility is that this is a unary percentage.
+      // If it succeeds, then we presume that this must be binary modulo, since the
+      // only things that parseUnary can handle are _higher_ precedence than unary %.
+      try {
+        parseUnary(state)
+        // Not sure if we could somehow use the result of that parseUnary? Without
+        // further analysis/testing, safer just to discard and let the parse proceed
+        Object.assign(state, previousState)
+      } catch {
+        // Not seeing a term at this point, so was a unary %
+        node = new OperatorNode('/', 'divide', [node, new ConstantNode(100)], false, true)
+      }
+    }
+
+    return node
+  }
+
   /**
    * Unary plus and minus, and logical and bitwise not
    * @return {Node} node
diff --git a/src/function/algebra/sylvester.js b/src/function/algebra/sylvester.js
index ab882374a8..53d0c3d0b5 100644
--- a/src/function/algebra/sylvester.js
+++ b/src/function/algebra/sylvester.js
@@ -35,7 +35,8 @@ export const createSylvester = /* #__PURE__ */ factory(name, dependencies, (
     subtract,
     identity,
     lusolve,
-    abs
+    abs,
+    config
   }
 ) => {
   /**
@@ -110,7 +111,7 @@ export const createSylvester = /* #__PURE__ */ factory(name, dependencies, (
 
     for (let k = 0; k < n; k++) {
       if (k < (n - 1) && abs(subset(G, index(k + 1, k))) > 1e-5) {
-        let RHS = vc(subset(D, index(all, k)), subset(D, index(all, k + 1)))
+        let RHS = vc(subset(D, index(all, [k])), subset(D, index(all, [k + 1])))
         for (let j = 0; j < k; j++) {
           RHS = add(RHS,
             vc(multiply(y[j], subset(G, index(j, k))), multiply(y[j], subset(G, index(j, k + 1))))
@@ -125,11 +126,11 @@ export const createSylvester = /* #__PURE__ */ factory(name, dependencies, (
           hc(gkm, add(F, gmm))
         )
         const yAux = lusolve(LHS, RHS)
-        y[k] = yAux.subset(index(range(0, m), 0))
-        y[k + 1] = yAux.subset(index(range(m, 2 * m), 0))
+        y[k] = yAux.subset(index(range(0, m), [0]))
+        y[k + 1] = yAux.subset(index(range(m, 2 * m), [0]))
         k++
       } else {
-        let RHS = subset(D, index(all, k))
+        let RHS = subset(D, index(all, [k]))
         for (let j = 0; j < k; j++) { RHS = add(RHS, multiply(y[j], subset(G, index(j, k)))) }
         const gkk = subset(G, index(k, k))
         const LHS = subtract(F, multiply(gkk, identity(m)))
@@ -139,7 +140,6 @@ export const createSylvester = /* #__PURE__ */ factory(name, dependencies, (
     }
     const Y = matrix(matrixFromColumns(...y))
     const X = multiply(U, multiply(Y, transpose(V)))
-
     return X
   }
 })
diff --git a/src/function/matrix/column.js b/src/function/matrix/column.js
index 971463729b..ef5e6c9dc3 100644
--- a/src/function/matrix/column.js
+++ b/src/function/matrix/column.js
@@ -51,8 +51,9 @@ export const createColumn = /* #__PURE__ */ factory(name, dependencies, ({ typed
     validateIndex(column, value.size()[1])
 
     const rowRange = range(0, value.size()[0])
-    const index = new Index(rowRange, column)
+    const index = new Index(rowRange, [column])
     const result = value.subset(index)
+    // once config.legacySubset just return result
     return isMatrix(result)
       ? result
       : matrix([[result]])
diff --git a/src/function/matrix/kron.js b/src/function/matrix/kron.js
index 1583672cb6..7378954b6b 100644
--- a/src/function/matrix/kron.js
+++ b/src/function/matrix/kron.js
@@ -22,7 +22,7 @@ export const createKron = /* #__PURE__ */ factory(name, dependencies, ({ typed,
      *    // returns [ [ 1, 2, 0, 0 ], [ 3, 4, 0, 0 ], [ 0, 0, 1, 2 ], [ 0, 0, 3, 4 ] ]
      *
      *    math.kron([1,1], [2,3,4])
-     *    // returns [ [ 2, 3, 4, 2, 3, 4 ] ]
+     *    // returns [2, 3, 4, 2, 3, 4]
      *
      * See also:
      *
@@ -49,39 +49,38 @@ export const createKron = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   })
 
   /**
-     * Calculate the Kronecker product of two matrices / vectors
-     * @param {Array} a  First vector
-     * @param {Array} b  Second vector
-     * @returns {Array} Returns the Kronecker product of x and y
-     * @private
+   * Calculate the Kronecker product of two (1-dimensional) vectors,
+   * with no dimension checking
+   * @param {Array} a  First vector
+   * @param {Array} b  Second vector
+   * @returns {Array}  the 1-dimensional Kronecker product of a and b
+   * @private
+   */
+  function _kron1d (a, b) {
+    // TODO in core overhaul: would be faster to see if we can choose a
+    // particular implementation of multiplyScalar at the beginning,
+    // rather than re-dispatch for _every_ ordered pair of entries.
+    return a.flatMap(x => b.map(y => multiplyScalar(x, y)))
+  }
+
+  /**
+   * Calculate the Kronecker product of two possibly multidimensional arrays
+   * @param {Array} a  First array
+   * @param {Array} b  Second array
+   * @param {number} [d]  common dimension; if missing, compute and match args
+   * @returns {Array} Returns the Kronecker product of x and y
+   * @private
      */
-  function _kron (a, b) {
-    // Deal with the dimensions of the matricies.
-    if (size(a).length === 1) {
-      // Wrap it in a 2D Matrix
-      a = [a]
-    }
-    if (size(b).length === 1) {
-      // Wrap it in a 2D Matrix
-      b = [b]
-    }
-    if (size(a).length > 2 || size(b).length > 2) {
-      throw new RangeError('Vectors with dimensions greater then 2 are not supported expected ' +
-            '(Size x = ' + JSON.stringify(a.length) + ', y = ' + JSON.stringify(b.length) + ')')
+  function _kron (a, b, d = -1) {
+    if (d < 0) {
+      let adim = size(a).length
+      let bdim = size(b).length
+      d = Math.max(adim, bdim)
+      while (adim++ < d) a = [a]
+      while (bdim++ < d) b = [b]
     }
-    const t = []
-    let r = []
 
-    return a.map(function (a) {
-      return b.map(function (b) {
-        r = []
-        t.push(r)
-        return a.map(function (y) {
-          return b.map(function (x) {
-            return r.push(multiplyScalar(y, x))
-          })
-        })
-      })
-    }) && t
+    if (d === 1) return _kron1d(a, b)
+    return a.flatMap(aSlice => b.map(bSlice => _kron(aSlice, bSlice, d - 1)))
   }
 })
diff --git a/src/function/matrix/row.js b/src/function/matrix/row.js
index e4a25738b5..5e99229404 100644
--- a/src/function/matrix/row.js
+++ b/src/function/matrix/row.js
@@ -51,8 +51,9 @@ export const createRow = /* #__PURE__ */ factory(name, dependencies, ({ typed, I
     validateIndex(row, value.size()[0])
 
     const columnRange = range(0, value.size()[1])
-    const index = new Index(row, columnRange)
+    const index = new Index([row], columnRange)
     const result = value.subset(index)
+    // once config.legacySubset just return result
     return isMatrix(result)
       ? result
       : matrix([[result]])
diff --git a/src/function/matrix/subset.js b/src/function/matrix/subset.js
index 48d7571315..8ca59e0523 100644
--- a/src/function/matrix/subset.js
+++ b/src/function/matrix/subset.js
@@ -21,8 +21,8 @@ export const createSubset = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *     // get a subset
    *     const d = [[1, 2], [3, 4]]
    *     math.subset(d, math.index(1, 0))             // returns 3
-   *     math.subset(d, math.index([0, 1], 1))        // returns [[2], [4]]
-   *     math.subset(d, math.index([false, true], 0)) // returns [[3]]
+   *     math.subset(d, math.index([0, 1], [1]))        // returns [[2], [4]]
+   *     math.subset(d, math.index([false, true], [0])) // returns [[3]]
    *
    *     // replace a subset
    *     const e = []
@@ -32,9 +32,9 @@ export const createSubset = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *
    *     // get submatrix using ranges
    *     const M = [
-   *       [1,2,3],
-   *       [4,5,6],
-   *       [7,8,9]
+   *       [1, 2, 3],
+   *       [4, 5, 6],
+   *       [7, 8, 9]
    *     ]
    *     math.subset(M, math.index(math.range(0,2), math.range(0,3))) // [[1, 2, 3], [4, 5, 6]]
    *
@@ -115,7 +115,7 @@ export const createSubset = /* #__PURE__ */ factory(name, dependencies, ({ typed
     if (typeof replacement === 'string') {
       throw new Error('can\'t boradcast a string')
     }
-    if (index._isScalar) {
+    if (index.isScalar()) {
       return replacement
     }
 
@@ -160,9 +160,14 @@ function _getSubstring (str, index) {
   const range = index.dimension(0)
 
   let substr = ''
-  range.forEach(function (v) {
+  function callback (v) {
     substr += str.charAt(v)
-  })
+  }
+  if (Number.isInteger(range)) {
+    callback(range)
+  } else {
+    range.forEach(callback)
+  }
 
   return substr
 }
@@ -196,7 +201,7 @@ function _setSubstring (str, index, replacement, defaultValue) {
   }
 
   const range = index.dimension(0)
-  const len = range.size()[0]
+  const len = Number.isInteger(range) ? 1 : range.size()[0]
 
   if (len !== replacement.length) {
     throw new DimensionError(range.size()[0], replacement.length)
@@ -213,9 +218,15 @@ function _setSubstring (str, index, replacement, defaultValue) {
     chars[i] = str.charAt(i)
   }
 
-  range.forEach(function (v, i) {
+  function callback (v, i) {
     chars[v] = replacement.charAt(i[0])
-  })
+  }
+
+  if (Number.isInteger(range)) {
+    callback(range, [0])
+  } else {
+    range.forEach(callback)
+  }
 
   // initialize undefined characters with a space
   if (chars.length > strLen) {
diff --git a/src/type/matrix/DenseMatrix.js b/src/type/matrix/DenseMatrix.js
index 9f1913ae77..112e0abc33 100644
--- a/src/type/matrix/DenseMatrix.js
+++ b/src/type/matrix/DenseMatrix.js
@@ -10,10 +10,11 @@ import { optimizeCallback } from '../../utils/optimizeCallback.js'
 
 const name = 'DenseMatrix'
 const dependencies = [
-  'Matrix'
+  'Matrix',
+  'config'
 ]
 
-export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies, ({ Matrix }) => {
+export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies, ({ Matrix, config }) => {
   /**
    * Dense Matrix implementation. A regular, dense matrix, supporting multi-dimensional matrices. This is the default matrix type.
    * @class DenseMatrix
@@ -218,7 +219,9 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
       throw new TypeError('Invalid index')
     }
 
-    const isScalar = index.isScalar()
+    const isScalar = config.legacySubset
+      ? index.size().every(idx => idx === 1)
+      : index.isScalar()
     if (isScalar) {
       // return a scalar
       return matrix.get(index.min())
@@ -238,12 +241,12 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
       }
 
       // retrieve submatrix
-      const returnMatrix = new DenseMatrix([])
+      const returnMatrix = new DenseMatrix()
       const submatrix = _getSubmatrix(matrix._data, index)
       returnMatrix._size = submatrix.size
       returnMatrix._datatype = matrix._datatype
       returnMatrix._data = submatrix.data
-      return returnMatrix
+      return config.legacySubset ? returnMatrix.reshape(index.size()) : returnMatrix
     }
   }
 
@@ -259,21 +262,31 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
   function _getSubmatrix (data, index) {
     const maxDepth = index.size().length - 1
     const size = Array(maxDepth)
-    return { data: getSubmatrixRecursive(data), size }
+    return { data: getSubmatrixRecursive(data), size: size.filter(x => x !== null) }
 
     function getSubmatrixRecursive (data, depth = 0) {
-      const ranges = index.dimension(depth)
-      size[depth] = ranges.size()[0]
+      const dims = index.dimension(depth)
+      function _mapIndex (dim, callback) {
+        // applies a callback for when the index is a Number or a Matrix
+        if (isNumber(dim)) return callback(dim)
+        else return dim.map(callback).valueOf()
+      }
+
+      if (isNumber(dims)) {
+        size[depth] = null
+      } else {
+        size[depth] = dims.size()[0]
+      }
       if (depth < maxDepth) {
-        return ranges.map(rangeIndex => {
-          validateIndex(rangeIndex, data.length)
-          return getSubmatrixRecursive(data[rangeIndex], depth + 1)
-        }).valueOf()
+        return _mapIndex(dims, dimIndex => {
+          validateIndex(dimIndex, data.length)
+          return getSubmatrixRecursive(data[dimIndex], depth + 1)
+        })
       } else {
-        return ranges.map(rangeIndex => {
-          validateIndex(rangeIndex, data.length)
-          return data[rangeIndex]
-        }).valueOf()
+        return _mapIndex(dims, dimIndex => {
+          validateIndex(dimIndex, data.length)
+          return data[dimIndex]
+        })
       }
     }
   }
@@ -300,19 +313,19 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
     const isScalar = index.isScalar()
 
     // calculate the size of the submatrix, and convert it into an Array if needed
-    let sSize
+    let submatrixSize
     if (isMatrix(submatrix)) {
-      sSize = submatrix.size()
+      submatrixSize = submatrix.size()
       submatrix = submatrix.valueOf()
     } else {
-      sSize = arraySize(submatrix)
+      submatrixSize = arraySize(submatrix)
     }
 
     if (isScalar) {
       // set a scalar
 
       // check whether submatrix is a scalar
-      if (sSize.length !== 0) {
+      if (submatrixSize.length !== 0) {
         throw new TypeError('Scalar expected')
       }
       matrix.set(index.min(), submatrix, defaultValue)
@@ -320,16 +333,16 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
       // set a submatrix
 
       // broadcast submatrix
-      if (!deepStrictEqual(sSize, iSize)) {
-        try {
-          if (sSize.length === 0) {
-            submatrix = broadcastTo([submatrix], iSize)
-          } else {
+      if (!deepStrictEqual(submatrixSize, iSize)) {
+        if (submatrixSize.length === 0) {
+          submatrix = broadcastTo([submatrix], iSize)
+        } else {
+          try {
             submatrix = broadcastTo(submatrix, iSize)
+          } catch (error) {
           }
-          sSize = arraySize(submatrix)
-        } catch {
         }
+        submatrixSize = arraySize(submatrix)
       }
 
       // validate dimensions
@@ -337,11 +350,11 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
         throw new DimensionError(iSize.length, matrix._size.length, '<')
       }
 
-      if (sSize.length < iSize.length) {
+      if (submatrixSize.length < iSize.length) {
         // calculate number of missing outer dimensions
         let i = 0
         let outer = 0
-        while (iSize[i] === 1 && sSize[i] === 1) {
+        while (iSize[i] === 1 && submatrixSize[i] === 1) {
           i++
         }
         while (iSize[i] === 1) {
@@ -350,12 +363,12 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
         }
 
         // unsqueeze both outer and inner dimensions
-        submatrix = unsqueeze(submatrix, iSize.length, outer, sSize)
+        submatrix = unsqueeze(submatrix, iSize.length, outer, submatrixSize)
       }
 
       // check whether the size of the submatrix matches the index size
-      if (!deepStrictEqual(iSize, sSize)) {
-        throw new DimensionError(iSize, sSize, '>')
+      if (!deepStrictEqual(iSize, submatrixSize)) {
+        throw new DimensionError(iSize, submatrixSize, '>')
       }
 
       // enlarge matrix when needed
@@ -386,16 +399,21 @@ export const createDenseMatrixClass = /* #__PURE__ */ factory(name, dependencies
 
     function setSubmatrixRecursive (data, submatrix, depth = 0) {
       const range = index.dimension(depth)
+      const recursiveCallback = (rangeIndex, i) => {
+        validateIndex(rangeIndex, data.length)
+        setSubmatrixRecursive(data[rangeIndex], submatrix[i[0]], depth + 1)
+      }
+      const finalCallback = (rangeIndex, i) => {
+        validateIndex(rangeIndex, data.length)
+        data[rangeIndex] = submatrix[i[0]]
+      }
+
       if (depth < maxDepth) {
-        range.forEach((rangeIndex, i) => {
-          validateIndex(rangeIndex, data.length)
-          setSubmatrixRecursive(data[rangeIndex], submatrix[i[0]], depth + 1)
-        })
+        if (isNumber(range)) recursiveCallback(range, [0])
+        else range.forEach(recursiveCallback)
       } else {
-        range.forEach((rangeIndex, i) => {
-          validateIndex(rangeIndex, data.length)
-          data[rangeIndex] = submatrix[i[0]]
-        })
+        if (isNumber(range)) finalCallback(range, [0])
+        else range.forEach(finalCallback)
       }
     }
   }
diff --git a/src/type/matrix/MatrixIndex.js b/src/type/matrix/MatrixIndex.js
index a76ab8526b..00f5447349 100644
--- a/src/type/matrix/MatrixIndex.js
+++ b/src/type/matrix/MatrixIndex.js
@@ -1,4 +1,4 @@
-import { isArray, isMatrix, isRange } from '../../utils/is.js'
+import { isArray, isMatrix, isRange, isNumber, isString } from '../../utils/is.js'
 import { clone } from '../../utils/object.js'
 import { isInteger } from '../../utils/number.js'
 import { factory } from '../../utils/factory.js'
@@ -29,7 +29,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
    * @Constructor Index
    * @param {...*} ranges
    */
-  function Index (ranges) {
+  function Index (...ranges) {
     if (!(this instanceof Index)) {
       throw new SyntaxError('Constructor must be called with the new operator')
     }
@@ -38,8 +38,8 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
     this._sourceSize = []
     this._isScalar = true
 
-    for (let i = 0, ii = arguments.length; i < ii; i++) {
-      const arg = arguments[i]
+    for (let i = 0, ii = ranges.length; i < ii; i++) {
+      const arg = ranges[i]
       const argIsArray = isArray(arg)
       const argIsMatrix = isMatrix(arg)
       const argType = typeof arg
@@ -50,6 +50,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
       } else if (argIsArray || argIsMatrix) {
         // create matrix
         let m
+        this._isScalar = false
 
         if (getMatrixDataType(arg) === 'boolean') {
           if (argIsArray) m = _createImmutableMatrix(_booleansArrayToNumbersForIndex(arg).valueOf())
@@ -60,16 +61,10 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
         }
 
         this._dimensions.push(m)
-        // size
-        const size = m.size()
-        // scalar
-        if (size.length !== 1 || size[0] !== 1 || sourceSize !== null) {
-          this._isScalar = false
-        }
       } else if (argType === 'number') {
-        this._dimensions.push(_createImmutableMatrix([arg]))
+        this._dimensions.push(arg)
       } else if (argType === 'bigint') {
-        this._dimensions.push(_createImmutableMatrix([Number(arg)]))
+        this._dimensions.push(Number(arg))
       } else if (argType === 'string') {
         // object property (arguments.count should be 1)
         this._dimensions.push(arg)
@@ -90,12 +85,15 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
   function _createImmutableMatrix (arg) {
     // loop array elements
     for (let i = 0, l = arg.length; i < l; i++) {
-      if (typeof arg[i] !== 'number' || !isInteger(arg[i])) {
+      if (!isNumber(arg[i]) || !isInteger(arg[i])) {
         throw new TypeError('Index parameters must be positive integer numbers')
       }
     }
     // create matrix
-    return new ImmutableDenseMatrix(arg)
+    const matrix = new ImmutableDenseMatrix()
+    matrix._data = arg
+    matrix._size = [arg.length]
+    return matrix
   }
 
   /**
@@ -134,7 +132,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
 
     for (let i = 0, ii = this._dimensions.length; i < ii; i++) {
       const d = this._dimensions[i]
-      size[i] = (typeof d === 'string') ? 1 : d.size()[0]
+      size[i] = (isString(d) || isNumber(d)) ? 1 : d.size()[0]
     }
 
     return size
@@ -150,7 +148,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
 
     for (let i = 0, ii = this._dimensions.length; i < ii; i++) {
       const range = this._dimensions[i]
-      values[i] = (typeof range === 'string') ? range : range.max()
+      values[i] = (isString(range) || isNumber(range)) ? range : range.max()
     }
 
     return values
@@ -166,7 +164,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
 
     for (let i = 0, ii = this._dimensions.length; i < ii; i++) {
       const range = this._dimensions[i]
-      values[i] = (typeof range === 'string') ? range : range.min()
+      values[i] = (isString(range) || isNumber(range)) ? range : range.min()
     }
 
     return values
@@ -192,11 +190,11 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
    * @returns {Range | null} range
    */
   Index.prototype.dimension = function (dim) {
-    if (typeof dim !== 'number') {
+    if (!isNumber(dim)) {
       return null
     }
 
-    return this._dimensions[dim] || null
+    return this._dimensions[dim] ?? null
   }
 
   /**
@@ -204,7 +202,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
    * @returns {boolean} Returns true if the index is an object property
    */
   Index.prototype.isObjectProperty = function () {
-    return this._dimensions.length === 1 && typeof this._dimensions[0] === 'string'
+    return this._dimensions.length === 1 && isString(this._dimensions[0])
   }
 
   /**
@@ -238,7 +236,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
     const array = []
     for (let i = 0, ii = this._dimensions.length; i < ii; i++) {
       const dimension = this._dimensions[i]
-      array.push((typeof dimension === 'string') ? dimension : dimension.toArray())
+      array.push(isString(dimension) || isNumber(dimension) ? dimension : dimension.toArray())
     }
     return array
   }
@@ -261,7 +259,7 @@ export const createIndexClass = /* #__PURE__ */ factory(name, dependencies, ({ I
 
     for (let i = 0, ii = this._dimensions.length; i < ii; i++) {
       const dimension = this._dimensions[i]
-      if (typeof dimension === 'string') {
+      if (isString(dimension)) {
         strings.push(JSON.stringify(dimension))
       } else {
         strings.push(dimension.toString())
diff --git a/src/type/matrix/SparseMatrix.js b/src/type/matrix/SparseMatrix.js
index b0301ecd5b..bbb6c80696 100644
--- a/src/type/matrix/SparseMatrix.js
+++ b/src/type/matrix/SparseMatrix.js
@@ -295,12 +295,14 @@ export const createSparseMatrixClass = /* #__PURE__ */ factory(name, dependencie
     const pv = []
 
     // loop rows in resulting matrix
-    rows.forEach(function (i, r) {
+    function rowsCallback (i, r) {
       // update permutation vector
       pv[i] = r[0]
       // mark i in workspace
       w[i] = true
-    })
+    }
+    if (Number.isInteger(rows)) rowsCallback(rows, [0])
+    else rows.forEach(rowsCallback)
 
     // result matrix arrays
     const values = mvalues ? [] : undefined
@@ -308,7 +310,7 @@ export const createSparseMatrixClass = /* #__PURE__ */ factory(name, dependencie
     const ptr = []
 
     // loop columns in result matrix
-    columns.forEach(function (j) {
+    function columnsCallback (j) {
       // update ptr
       ptr.push(index.length)
       // loop values in column j
@@ -323,7 +325,9 @@ export const createSparseMatrixClass = /* #__PURE__ */ factory(name, dependencie
           if (values) { values.push(mvalues[k]) }
         }
       }
-    })
+    }
+    if (Number.isInteger(columns)) columnsCallback(columns)
+    else columns.forEach(columnsCallback)
     // update ptr
     ptr.push(index.length)
 
@@ -398,7 +402,7 @@ export const createSparseMatrixClass = /* #__PURE__ */ factory(name, dependencie
       if (iSize.length === 1) {
         // if the replacement index only has 1 dimension, go trough each one and set its value
         const range = index.dimension(0)
-        range.forEach(function (dataIndex, subIndex) {
+        _forEachIndex(range, (dataIndex, subIndex) => {
           validateIndex(dataIndex)
           matrix.set([dataIndex, 0], submatrix[subIndex[0]], defaultValue)
         })
@@ -406,9 +410,9 @@ export const createSparseMatrixClass = /* #__PURE__ */ factory(name, dependencie
         // if the replacement index has 2 dimensions, go through each one and set the value in the correct index
         const firstDimensionRange = index.dimension(0)
         const secondDimensionRange = index.dimension(1)
-        firstDimensionRange.forEach(function (firstDataIndex, firstSubIndex) {
+        _forEachIndex(firstDimensionRange, (firstDataIndex, firstSubIndex) => {
           validateIndex(firstDataIndex)
-          secondDimensionRange.forEach(function (secondDataIndex, secondSubIndex) {
+          _forEachIndex(secondDimensionRange, (secondDataIndex, secondSubIndex) => {
             validateIndex(secondDataIndex)
             matrix.set([firstDataIndex, secondDataIndex], submatrix[firstSubIndex[0]][secondSubIndex[0]], defaultValue)
           })
@@ -416,6 +420,12 @@ export const createSparseMatrixClass = /* #__PURE__ */ factory(name, dependencie
       }
     }
     return matrix
+
+    function _forEachIndex (index, callback) {
+      // iterate cases where index is a Matrix or a Number
+      if (isNumber(index)) callback(index, [0])
+      else index.forEach(callback)
+    }
   }
 
   /**
diff --git a/test/node-tests/defaultInstance.test.cjs b/test/node-tests/defaultInstance.test.cjs
index 350fb65f02..0a8eb7a41a 100644
--- a/test/node-tests/defaultInstance.test.cjs
+++ b/test/node-tests/defaultInstance.test.cjs
@@ -15,6 +15,7 @@ describe('defaultInstance', function () {
       predictable: false,
       relTol: 1e-12,
       absTol: 1e-15,
+      legacySubset: false,
       randomSeed: null
     })
   })
@@ -34,6 +35,7 @@ describe('defaultInstance', function () {
       predictable: false,
       relTol: 1e-12,
       absTol: 1e-15,
+      legacySubset: false,
       randomSeed: null
     })
   })
@@ -79,6 +81,7 @@ describe('defaultInstance', function () {
       predictable: true,
       relTol: 1e-12,
       absTol: 1e-15,
+      legacySubset: false,
       randomSeed: null
     })
 
@@ -101,6 +104,7 @@ describe('defaultInstance', function () {
       predictable: false,
       relTol: 1e-12,
       absTol: 1e-15,
+      legacySubset: false,
       randomSeed: null
     })
 
diff --git a/test/unit-tests/core/config.test.js b/test/unit-tests/core/config.test.js
index ac8b1aaa98..bfe730ced1 100644
--- a/test/unit-tests/core/config.test.js
+++ b/test/unit-tests/core/config.test.js
@@ -37,4 +37,28 @@ describe('config', function () {
     // Restore console.warn
     warnStub.restore()
   })
+
+  it('should work with config legacySubset during deprecation', function () {
+    const math2 = math.create()
+    // Add a spy to temporarily disable console.warn
+    const warnStub = sinon.stub(console, 'warn')
+
+    // Set legacySubset to true and should throw a warning
+    assert.doesNotThrow(function () { math2.config({ legacySubset: true }) })
+
+    // Check if legacySubset is set
+    assert.strictEqual(math2.config().legacySubset, true)
+
+    // Check if console.warn was called
+    assert.strictEqual(warnStub.callCount, 1)
+
+    // Set legacySubset to false, should not throw a warning
+    assert.doesNotThrow(function () { math2.config({ legacySubset: false }) })
+
+    // Validate that  if console.warn was not called again
+    assert.strictEqual(warnStub.callCount, 1)
+
+    // Restore console.warn
+    warnStub.restore()
+  })
 })
diff --git a/test/unit-tests/expression/node/AccessorNode.test.js b/test/unit-tests/expression/node/AccessorNode.test.js
index 675e362816..292361ee1d 100644
--- a/test/unit-tests/expression/node/AccessorNode.test.js
+++ b/test/unit-tests/expression/node/AccessorNode.test.js
@@ -73,7 +73,7 @@ describe('AccessorNode', function () {
     const scope = {
       a: [[1, 2], [3, 4]]
     }
-    assert.deepStrictEqual(expr.evaluate(scope), [[3, 4]])
+    assert.deepStrictEqual(expr.evaluate(scope), [3, 4])
   })
 
   it('should compile a AccessorNode with "end" in an expression', function () {
@@ -185,7 +185,7 @@ describe('AccessorNode', function () {
     const scope = {
       a: [[1, 2], [3, 4]]
     }
-    assert.deepStrictEqual(expr.evaluate(scope), [[4, 3]])
+    assert.deepStrictEqual(expr.evaluate(scope), [4, 3])
   })
 
   it('should compile a AccessorNode with "end" both as value and in a range', function () {
@@ -203,7 +203,7 @@ describe('AccessorNode', function () {
     const scope = {
       a: [[1, 2], [3, 4]]
     }
-    assert.deepStrictEqual(expr.evaluate(scope), [[3, 4]])
+    assert.deepStrictEqual(expr.evaluate(scope), [3, 4])
   })
 
   it('should use the inner context when using "end" in a nested index', function () {
diff --git a/test/unit-tests/expression/parse.test.js b/test/unit-tests/expression/parse.test.js
index 09455c5dc0..2da13b0552 100644
--- a/test/unit-tests/expression/parse.test.js
+++ b/test/unit-tests/expression/parse.test.js
@@ -289,6 +289,21 @@ describe('parse', function () {
       assert.strictEqual(parseAndEval('0x1.'), 1)
     })
 
+    it('should require hex, bin, oct values to be followed by whitespace or a delimiter', function () {
+      assert.throws(() => parseAndEval('0b0a'), /SyntaxError: String "0b0a" is not a valid number/)
+      assert.throws(() => parseAndEval('0x1k'), /SyntaxError: String "0x1k" is not a valid number/)
+      assert.throws(() => parseAndEval('0o1k'), /SyntaxError: String "0o1k" is not a valid number/)
+      assert.throws(() => parseAndEval('0b1k'), /SyntaxError: String "0b1k" is not a valid number/)
+
+      assert.strictEqual(parseAndEval('0x1 k', { k: 2 }), 2)
+      assert.strictEqual(parseAndEval('0o1 k', { k: 2 }), 2)
+      assert.strictEqual(parseAndEval('0b1 k', { k: 2 }), 2)
+
+      assert.strictEqual(parseAndEval('0x1*k', { k: 2 }), 2)
+      assert.strictEqual(parseAndEval('0o1*k', { k: 2 }), 2)
+      assert.strictEqual(parseAndEval('0b1*k', { k: 2 }), 2)
+    })
+
     it('should parse a number followed by e', function () {
       approxEqual(parseAndEval('2e'), 2 * Math.E)
     })
@@ -337,7 +352,7 @@ describe('parse', function () {
 
       assert.throws(function () { parseAndEval('0b123.45') }, /SyntaxError: String "0b123\.45" is not a valid number/)
       assert.throws(function () { parseAndEval('0o89.89') }, /SyntaxError: String "0o89\.89" is not a valid number/)
-      assert.throws(function () { parseAndEval('0xghji.xyz') }, /SyntaxError: String "0x" is not a valid number/)
+      assert.throws(function () { parseAndEval('0xghji.xyz') }, /SyntaxError: String "0xghji.xyz" is not a valid number/)
     })
   })
 
@@ -700,18 +715,18 @@ describe('parse', function () {
           [7, 8, 9]
         ])
       }
-      assert.deepStrictEqual(parseAndEval('a[2, :]', scope), math.matrix([[4, 5, 6]]))
-      assert.deepStrictEqual(parseAndEval('a[2, :2]', scope), math.matrix([[4, 5]]))
-      assert.deepStrictEqual(parseAndEval('a[2, :end-1]', scope), math.matrix([[4, 5]]))
-      assert.deepStrictEqual(parseAndEval('a[2, 2:]', scope), math.matrix([[5, 6]]))
-      assert.deepStrictEqual(parseAndEval('a[2, 2:3]', scope), math.matrix([[5, 6]]))
-      assert.deepStrictEqual(parseAndEval('a[2, 1:2:3]', scope), math.matrix([[4, 6]]))
-      assert.deepStrictEqual(parseAndEval('a[:, 2]', scope), math.matrix([[2], [5], [8]]))
-      assert.deepStrictEqual(parseAndEval('a[:2, 2]', scope), math.matrix([[2], [5]]))
-      assert.deepStrictEqual(parseAndEval('a[:end-1, 2]', scope), math.matrix([[2], [5]]))
-      assert.deepStrictEqual(parseAndEval('a[2:, 2]', scope), math.matrix([[5], [8]]))
-      assert.deepStrictEqual(parseAndEval('a[2:3, 2]', scope), math.matrix([[5], [8]]))
-      assert.deepStrictEqual(parseAndEval('a[1:2:3, 2]', scope), math.matrix([[2], [8]]))
+      assert.deepStrictEqual(parseAndEval('a[2, :]', scope), math.matrix([4, 5, 6]))
+      assert.deepStrictEqual(parseAndEval('a[2, :2]', scope), math.matrix([4, 5]))
+      assert.deepStrictEqual(parseAndEval('a[2, :end-1]', scope), math.matrix([4, 5]))
+      assert.deepStrictEqual(parseAndEval('a[2, 2:]', scope), math.matrix([5, 6]))
+      assert.deepStrictEqual(parseAndEval('a[2, 2:3]', scope), math.matrix([5, 6]))
+      assert.deepStrictEqual(parseAndEval('a[2, 1:2:3]', scope), math.matrix([4, 6]))
+      assert.deepStrictEqual(parseAndEval('a[:, 2]', scope), math.matrix([2, 5, 8]))
+      assert.deepStrictEqual(parseAndEval('a[:2, [2]]', scope), math.matrix([[2], [5]]))
+      assert.deepStrictEqual(parseAndEval('a[:end-1, [2]]', scope), math.matrix([[2], [5]]))
+      assert.deepStrictEqual(parseAndEval('a[2:, [2]]', scope), math.matrix([[5], [8]]))
+      assert.deepStrictEqual(parseAndEval('a[2:3, [2]]', scope), math.matrix([[5], [8]]))
+      assert.deepStrictEqual(parseAndEval('a[1:2:3, [2]]', scope), math.matrix([[2], [8]]))
     })
 
     it('should get a matrix subset of a matrix subset', function () {
@@ -722,7 +737,7 @@ describe('parse', function () {
           [7, 8, 9]
         ])
       }
-      assert.deepStrictEqual(parseAndEval('a[2, :][1,1]', scope), 4)
+      assert.deepStrictEqual(parseAndEval('a[[2], :][1,1]', scope), 4)
     })
 
     it('should get BigNumber value from an array', function () {
@@ -771,7 +786,7 @@ describe('parse', function () {
       assert.deepStrictEqual(parseAndEval('a[1:3,1:2]', scope), math.matrix([[100, 2], [3, 10], [0, 12]]))
 
       scope.b = [[1, 2], [3, 4]]
-      assert.deepStrictEqual(parseAndEval('b[1,:]', scope), [[1, 2]])
+      assert.deepStrictEqual(parseAndEval('b[1,:]', scope), [1, 2])
     })
 
     it('should get/set the matrix correctly for 3d matrices', function () {
@@ -792,10 +807,10 @@ describe('parse', function () {
       ]))
 
       assert.deepStrictEqual(parseAndEval('size(f)', scope), math.matrix([2, 2, 2], 'dense', 'number'))
-      assert.deepStrictEqual(parseAndEval('f[:,:,1]', scope), math.matrix([[[1], [2]], [[3], [4]]]))
-      assert.deepStrictEqual(parseAndEval('f[:,:,2]', scope), math.matrix([[[5], [6]], [[7], [8]]]))
-      assert.deepStrictEqual(parseAndEval('f[:,2,:]', scope), math.matrix([[[2, 6]], [[4, 8]]]))
-      assert.deepStrictEqual(parseAndEval('f[2,:,:]', scope), math.matrix([[[3, 7], [4, 8]]]))
+      assert.deepStrictEqual(parseAndEval('f[:,:,1]', scope), math.matrix([[1, 2], [3, 4]]))
+      assert.deepStrictEqual(parseAndEval('f[:,:,2]', scope), math.matrix([[5, 6], [7, 8]]))
+      assert.deepStrictEqual(parseAndEval('f[:,2,:]', scope), math.matrix([[2, 6], [4, 8]]))
+      assert.deepStrictEqual(parseAndEval('f[2,:,:]', scope), math.matrix([[3, 7], [4, 8]]))
 
       parseAndEval('a=diag([1,2,3,4])', scope)
       assert.deepStrictEqual(parseAndEval('a[3:end, 3:end]', scope), math.matrix([[3, 0], [0, 4]]))
@@ -1374,36 +1389,62 @@ describe('parse', function () {
     it('should parse % with multiplication', function () {
       approxEqual(parseAndEval('100*50%'), 50)
       approxEqual(parseAndEval('50%*100'), 50)
-      assert.throws(function () { parseAndEval('50%(*100)') }, /Value expected/)
+      assert.throws(function () { parseAndEval('50%(*100)') }, SyntaxError)
     })
 
     it('should parse % with division', function () {
-      approxEqual(parseAndEval('100/50%'), 0.02) // should be treated as ((100/50)%)
-      approxEqual(parseAndEval('100/50%*2'), 0.04) // should be treated as ((100/50)%))×2
+      approxEqual(parseAndEval('100/50%'), 200) // should be treated as 100/(50%)
+      approxEqual(parseAndEval('100/50%*2'), 400) // should be treated as (100/(50%))×2
       approxEqual(parseAndEval('50%/100'), 0.005)
-      assert.throws(function () { parseAndEval('50%(/100)') }, /Value expected/)
+      approxEqual(parseAndEval('50%(13)'), 11) // should be treated as 50 % (13)
+      assert.throws(function () { parseAndEval('50%(/100)') }, SyntaxError)
+    })
+
+    it('should parse unary % before division, binary % with division', function () {
+      approxEqual(parseAndEval('10/200%%3'), 2) // should be treated as (10/(200%))%3
+    })
+
+    it('should reject repeated unary percentage operators', function () {
+      assert.throws(function () { math.parse('17%%') }, SyntaxError)
+      assert.throws(function () { math.parse('17%%*5') }, SyntaxError)
+      assert.throws(function () { math.parse('10/200%%%3') }, SyntaxError)
+    })
+
+    it('should parse unary % before division, binary % with division', function () {
+      approxEqual(parseAndEval('10/200%%3'), 2) // should be treated as (10/(200%))%3
     })
 
-    it('should parse % with addition', function () {
+    it('should reject repeated unary percentage operators', function () {
+      assert.throws(function () { math.parse('17%%') }, /Unexpected end of expression/)
+      assert.throws(function () { math.parse('17%%*5') }, /Value expected \(char 5\)/)
+      assert.throws(function () { math.parse('10/200%%%3') }, /Value expected \(char 9\)/)
+    })
+
+    it('should parse unary % with addition', function () {
       approxEqual(parseAndEval('100+3%'), 103)
-      approxEqual(parseAndEval('3%+100'), 100.03)
+      assert.strictEqual(parseAndEval('3%+100'), 3) // treat as 3 mod 100
     })
 
-    it('should parse % with subtraction', function () {
+    it('should parse unary % with subtraction', function () {
       approxEqual(parseAndEval('100-3%'), 97)
-      approxEqual(parseAndEval('3%-100'), -99.97)
+      assert.strictEqual(parseAndEval('3%-100'), -97) // treat as 3 mod -100
+    })
+
+    it('should parse binary % with bitwise negation', function () {
+      assert.strictEqual(parseAndEval('11%~1'), -1) // equivalent to 11 mod -2
+      assert.strictEqual(parseAndEval('11%~-3'), 1) // equivalent to 11 mod 2
     })
 
     it('should parse operator mod', function () {
       approxEqual(parseAndEval('8 mod 3'), 2)
     })
 
-    it('should give equal precedence to % and * operators', function () {
+    it('should give equal precedence to binary % and * operators', function () {
       approxEqual(parseAndStringifyWithParens('10 % 3 * 2'), '(10 % 3) * 2')
       approxEqual(parseAndStringifyWithParens('10 * 3 % 4'), '(10 * 3) % 4')
     })
 
-    it('should give equal precedence to % and / operators', function () {
+    it('should give equal precedence to binary % and / operators', function () {
       approxEqual(parseAndStringifyWithParens('10 % 4 / 2'), '(10 % 4) / 2')
       approxEqual(parseAndStringifyWithParens('10 / 2 % 3'), '(10 / 2) % 3')
     })
@@ -1418,12 +1459,12 @@ describe('parse', function () {
       approxEqual(parseAndStringifyWithParens('8 / 3 mod 2'), '(8 / 3) mod 2')
     })
 
-    it('should give equal precedence to % and .* operators', function () {
+    it('should give equal precedence to binary % and .* operators', function () {
       approxEqual(parseAndStringifyWithParens('10 % 3 .* 2'), '(10 % 3) .* 2')
       approxEqual(parseAndStringifyWithParens('10 .* 3 % 4'), '(10 .* 3) % 4')
     })
 
-    it('should give equal precedence to % and ./ operators', function () {
+    it('should give equal precedence to binary % and ./ operators', function () {
       approxEqual(parseAndStringifyWithParens('10 % 4 ./ 2'), '(10 % 4) ./ 2')
       approxEqual(parseAndStringifyWithParens('10 ./ 2 % 3'), '(10 ./ 2) % 3')
     })
diff --git a/test/unit-tests/function/algebra/sylvester.test.js b/test/unit-tests/function/algebra/sylvester.test.js
index 3af2472f9d..72bb42aa55 100644
--- a/test/unit-tests/function/algebra/sylvester.test.js
+++ b/test/unit-tests/function/algebra/sylvester.test.js
@@ -1,7 +1,7 @@
 // test schur decomposition
 import assert from 'assert'
-
 import math from '../../../../src/defaultInstance.js'
+import sinon from 'sinon'
 
 describe('sylvester', function () {
   it('should solve sylvester equation of order 5 with Arrays', function () {
@@ -53,4 +53,51 @@ describe('sylvester', function () {
       [-1.0935231387905004, 4.113817086842746, 5.747671819196675, -0.9408309030864932, 2.967655969930743]
     ])), 'fro') < 1e-3)
   })
+
+  it('should work with config legacySubset during deprecation', function () {
+    const math2 = math.create()
+    // Add a spy to temporarily disable console.warn
+    const warnStub = sinon.stub(console, 'warn')
+
+    math2.config({ legacySubset: true })
+
+    // Test legacy syntax with sylvester
+    // This is not strictly necessary and shoudl be removed after the deprecation period
+    const sylvesterA = [[-5.3, -1.4, -0.2, 0.7],
+      [-0.4, -1.0, -0.1, -1.2],
+      [0.3, 0.7, -2.5, 0.7],
+      [3.6, -0.1, 1.4, -2.4]]
+
+    const sylvesterB = [
+      [1.1, -0.3, -0.9, 0.8, -2.5],
+      [-0.6, 2.6, 0.2, 0.4, 1.3],
+      [0.4, -0.5, -0.2, 0.2, -0.1],
+      [-0.4, -1.9, -0.2, 0.5, 1.4],
+      [0.4, -1.0, -0.1, -0.8, -1.3]
+    ]
+    const sylvesterC = [
+      [1.4, 1.1, -1.9, 0.1, 1.2],
+      [-1.7, 0.1, -0.4, 2.1, 0.5],
+      [1.9, 2.3, -0.8, -0.7, 1.0],
+      [-1.1, 2.8, -0.8, -0.3, 0.9]
+    ]
+
+    const sylvesterSol = [
+      [-0.19393862606643053, -0.17101629636521865, 2.709348263225366, -0.0963000767188319, 0.5244718194343121],
+      [0.38421326955977486, -0.21159588555260944, -6.544262021555474, -0.15113424769761136, -2.312533293658291],
+      [-2.2708235174374747, 4.498279916441834, 1.4553799673144823, -0.9300926971755248, 2.5508111398452353],
+      [-1.0935231387905004, 4.113817086842746, 5.747671819196675, -0.9408309030864932, 2.967655969930743]
+    ]
+    assert.ok(math2.norm(math2.subtract(
+      math2.sylvester(sylvesterA, sylvesterB, sylvesterC), sylvesterSol), 'fro') < 1e-3)
+
+    // Test without legacySubset
+    math2.config({ legacySubset: false })
+
+    assert.ok(math2.norm(math2.subtract(
+      math2.sylvester(sylvesterA, sylvesterB, sylvesterC), sylvesterSol), 'fro') < 1e-3)
+
+    // Restore console.warn
+    warnStub.restore()
+  })
 })
diff --git a/test/unit-tests/function/matrix/column.test.js b/test/unit-tests/function/matrix/column.test.js
index d1e11d5725..6c6e3fed57 100644
--- a/test/unit-tests/function/matrix/column.test.js
+++ b/test/unit-tests/function/matrix/column.test.js
@@ -1,5 +1,6 @@
 import assert from 'assert'
 import math from '../../../../src/defaultInstance.js'
+import sinon from 'sinon'
 
 const column = math.column
 const matrix = math.matrix
@@ -107,4 +108,37 @@ describe('column', function () {
       c.valueOf(), [[0], [0], [0], [0], [0]]
     )
   })
+
+  it('should work with config legacySubset during deprecation', function () {
+    const math2 = math.create()
+    // Add a spy to temporarily disable console.warn
+    const warnStub = sinon.stub(console, 'warn')
+
+    math2.config({ legacySubset: true })
+
+    const a = [
+      [0, 2, 0, 0, 0],
+      [0, 1, 0, 2, 4],
+      [0, 0, 0, 0, 0],
+      [8, 4, 0, 3, 0],
+      [0, 0, 0, 6, 0]
+    ]
+    // Test column with legacySubset syntax
+    // This is not strictly necessary and shoudl be removed after the deprecation period
+
+    assert.deepStrictEqual(
+      math2.column(a, 4).valueOf(), [[0], [4], [0], [0], [0]]
+    )
+
+    // Test column with legacySubset syntax
+    math2.config({ legacySubset: false })
+
+    // Test column without legacySubset syntax
+    assert.deepStrictEqual(
+      math2.column(a, 4).valueOf(), [[0], [4], [0], [0], [0]]
+    )
+
+    // Restore console.warn
+    warnStub.restore()
+  })
 })
diff --git a/test/unit-tests/function/matrix/kron.test.js b/test/unit-tests/function/matrix/kron.test.js
index f93516c696..a344c52e62 100644
--- a/test/unit-tests/function/matrix/kron.test.js
+++ b/test/unit-tests/function/matrix/kron.test.js
@@ -5,6 +5,7 @@ import math from '../../../../src/defaultInstance.js'
 
 describe('kron', function () {
   it('should calculate the Kronecker product of two arrays', function () {
+    assert.deepStrictEqual(math.kron([[2]], [[3]]), [[6]])
     assert.deepStrictEqual(math.kron([
       [1, -2, 1],
       [1, 1, 0]
@@ -28,14 +29,27 @@ describe('kron', function () {
     ])
   })
 
+  it('should calculate product for empty 1D Arrays', function () {
+    assert.deepStrictEqual(math.kron([], []), [])
+  })
+
   it('should calculate product for empty 2D Arrays', function () {
     assert.deepStrictEqual(math.kron([[]], [[]]), [[]])
   })
 
   it('should calculate product for 1D Arrays', function () {
+    assert.deepStrictEqual(math.kron([2], [3]), [6])
+    assert.deepStrictEqual(math.kron([1, 2], [3, 4]), [3, 4, 6, 8])
+    assert.deepStrictEqual(math.kron([1, 2, 6, 8], [12, 1, 2, 3]), [12, 1, 2, 3, 24, 2, 4, 6, 72, 6, 12, 18, 96, 8, 16, 24])
+  })
+
+  it('should calculate product for 1D & 2D Arrays', function () {
     assert.deepStrictEqual(math.kron([1, 1], [[1, 0], [0, 1]]), [[1, 0, 1, 0], [0, 1, 0, 1]])
     assert.deepStrictEqual(math.kron([[1, 0], [0, 1]], [1, 1]), [[1, 1, 0, 0], [0, 0, 1, 1]])
-    assert.deepStrictEqual(math.kron([1, 2, 6, 8], [12, 1, 2, 3]), [[12, 1, 2, 3, 24, 2, 4, 6, 72, 6, 12, 18, 96, 8, 16, 24]])
+    assert.deepStrictEqual(math.kron([[1, 2]], [[1, 2, 3]]), [[1, 2, 3, 2, 4, 6]])
+    assert.deepStrictEqual(math.kron([[1], [2]], [[1], [2], [3]]), [[1], [2], [3], [2], [4], [6]])
+    assert.deepStrictEqual(math.kron([[1, 2]], [[1], [2], [3]]), [[1, 2], [2, 4], [3, 6]])
+    assert.deepStrictEqual(math.kron([[1], [2]], [[1, 2, 3]]), [[1, 2, 3], [2, 4, 6]])
   })
 
   it('should support complex numbers', function () {
@@ -55,10 +69,32 @@ describe('kron', function () {
     ])
   })
 
-  it('should throw an error for greater then 2 dimensions', function () {
-    assert.throws(function () {
-      math.kron([[[1, 1], [1, 1]], [[1, 1], [1, 1]]], [[[1, 2, 3], [4, 5, 6]], [[6, 7, 8], [9, 10, 11]]])
-    })
+  it('should calculate a 3D Kronecker product', function () {
+    assert.deepStrictEqual(
+      math.kron([
+        [[1, 2], [2, 3]],
+        [[2, 3], [3, 4]]
+      ], [
+        [[4, 3], [3, 2]],
+        [[3, 2], [2, 1]]
+      ]), [
+        /* eslint-disable no-multi-spaces, array-bracket-spacing */
+        [[4, 3,  8, 6], [3, 2, 6, 4], [ 8, 6, 12,  9], [6, 4,  9, 6]],
+        [[3, 2,  6, 4], [2, 1, 4, 2], [ 6, 4,  9,  6], [4, 2,  6, 3]],
+        [[8, 6, 12, 9], [6, 4, 9, 6], [12, 9, 16, 12], [9, 6, 12, 8]],
+        [[6, 4,  9, 6], [4, 2, 6, 3], [ 9, 6, 12,  8], [6, 3,  8, 4]]
+        /* eslint-enable */
+      ]
+    )
+  })
+
+  it('should allow mixed-dimensional Kronecker products', function () {
+    const b = [[[4, 3], [3, 2]], [[3, 2], [2, 1]]]
+    const a = [1, 2]
+    assert.deepStrictEqual(math.kron(a, b), math.kron([[a]], b))
+    assert.deepStrictEqual(math.kron([a], b), math.kron([[a]], b))
+    assert.deepStrictEqual(math.kron(b, a), math.kron(b, [[a]]))
+    assert.deepStrictEqual(math.kron(b, [a]), math.kron(b, [[a]]))
   })
 
   it('should throw an error if called with an invalid number of arguments', function () {
@@ -81,10 +117,12 @@ describe('kron', function () {
       assert.deepStrictEqual(product.toArray(), [[13, 26, 0, 0], [715, -13, 0, -0], [0, 0, -1, -2], [0, -0, -55, 1]])
     })
 
-    it('should throw an error for invalid Kronecker product of matrix', function () {
-      const y = math.matrix([[[]]])
+    it('should calculate the Kronecker product of 3d matrices', function () {
+      const y = math.matrix([[[3]]])
       const x = math.matrix([[[1, 1], [1, 1]], [[1, 1], [1, 1]]])
-      assert.throws(function () { math.kron(y, x) })
+      const product = math.kron(x, y)
+      assert.deepStrictEqual(
+        product.toArray(), [[[3, 3], [3, 3]], [[3, 3], [3, 3]]])
     })
   })
 
diff --git a/test/unit-tests/function/matrix/map.test.js b/test/unit-tests/function/matrix/map.test.js
index 2211f733ac..387f4e5f93 100644
--- a/test/unit-tests/function/matrix/map.test.js
+++ b/test/unit-tests/function/matrix/map.test.js
@@ -281,7 +281,7 @@ describe('map', function () {
   it('should operate from the parser with multiple inputs that need broadcasting and one based indices and the broadcasted arrays', function () {
     // this is a convoluted way of calculating f(a,b,idx) = 2a+2b+index
     // 2(1) + 2([3,4]) + [1, 2] # yields [9, 12]
-    const arr2 = math.evaluate('map([1],[3,4], f(a,b,idx,A,B)= a + A[idx] + b + B[idx] + idx[1])')
+    const arr2 = math.evaluate('map([1],[3,4], f(a,b,idx,A,B)= a + A[idx[1]] + b + B[idx[1]] + idx[1])')
     const expected = math.matrix([9, 12])
     assert.deepStrictEqual(arr2, expected)
   })
diff --git a/test/unit-tests/function/matrix/row.test.js b/test/unit-tests/function/matrix/row.test.js
index c864981ef0..1b0d23aa15 100644
--- a/test/unit-tests/function/matrix/row.test.js
+++ b/test/unit-tests/function/matrix/row.test.js
@@ -1,5 +1,6 @@
 import assert from 'assert'
 import math from '../../../../src/defaultInstance.js'
+import sinon from 'sinon'
 
 const row = math.row
 const matrix = math.matrix
@@ -107,4 +108,36 @@ describe('row', function () {
       r.valueOf(), [[0, 0, 0, 0, 0]]
     )
   })
+
+  it('should work with config legacySubset during deprecation', function () {
+    const math2 = math.create()
+    // Add a spy to temporarily disable console.warn
+    const warnStub = sinon.stub(console, 'warn')
+
+    math2.config({ legacySubset: true })
+
+    const a = [
+      [0, 2, 0, 0, 0],
+      [0, 1, 0, 2, 4],
+      [0, 0, 0, 0, 0],
+      [8, 4, 0, 3, 0],
+      [0, 0, 0, 6, 0]
+    ]
+
+    // Test row with legacySubset syntax
+    // This is not strictly necessary and should be removed after the deprecation period
+    assert.deepStrictEqual(
+      math2.row(a, 3).valueOf(), [[8, 4, 0, 3, 0]]
+    )
+
+    // Test row without legacySubset syntax
+    math2.config({ legacySubset: false })
+
+    assert.deepStrictEqual(
+      math2.row(a, 3).valueOf(), [[8, 4, 0, 3, 0]]
+    )
+
+    // Restore console.warn
+    warnStub.restore()
+  })
 })
diff --git a/test/unit-tests/function/matrix/subset.test.js b/test/unit-tests/function/matrix/subset.test.js
index 10b5003a1e..a99fcfcac4 100644
--- a/test/unit-tests/function/matrix/subset.test.js
+++ b/test/unit-tests/function/matrix/subset.test.js
@@ -1,6 +1,8 @@
 import assert from 'assert'
 import math from '../../../../src/defaultInstance.js'
 import { DimensionError } from '../../../../src/error/DimensionError.js'
+import sinon from 'sinon'
+
 const subset = math.subset
 const matrix = math.matrix
 const Range = math.Range
@@ -11,13 +13,13 @@ describe('subset', function () {
   const b = math.matrix(a)
 
   it('should get the right subset of an array', function () {
-    assert.deepStrictEqual(subset(a, index(new Range(0, 2), 1)), [[2], [4]])
+    assert.deepStrictEqual(subset(a, index(new Range(0, 2), 1)), [2, 4])
     assert.deepStrictEqual(subset(a, index(1, 0)), 3)
     assert.deepStrictEqual(subset([math.bignumber(2)], index(0)), math.bignumber(2))
   })
 
   it('should get the right subset of an array of booleans', function () {
-    assert.deepStrictEqual(subset(a, index([true, true], 1)), [[2], [4]])
+    assert.deepStrictEqual(subset(a, index([true, true], [1])), [[2], [4]])
     assert.deepStrictEqual(subset(a, index([false, true], [true, false])), [[3]])
     assert.deepStrictEqual(subset([math.bignumber(2)], index([true])), [math.bignumber(2)])
   })
@@ -32,7 +34,7 @@ describe('subset', function () {
   })
 
   it('should get the right subset of an array of booleans in the parser', function () {
-    assert.deepStrictEqual(math.evaluate('a[[true, true], 2]', { a }), [[2], [4]])
+    assert.deepStrictEqual(math.evaluate('a[[true, true], 2]', { a }), [2, 4])
     assert.deepStrictEqual(math.evaluate('a[[false, true], [true, false]]', { a }), [[3]])
     assert.deepStrictEqual(math.evaluate('[bignumber(2)][[true]]'), math.matrix([math.bignumber(2)]))
   })
@@ -75,7 +77,7 @@ describe('subset', function () {
   })
 
   it('should get the right subset of a matrix', function () {
-    assert.deepStrictEqual(subset(b, index(new Range(0, 2), 1)), matrix([[2], [4]]))
+    assert.deepStrictEqual(subset(b, index(new Range(0, 2), 1)), matrix([2, 4]))
     assert.deepStrictEqual(subset(b, index(1, 0)), 3)
   })
 
@@ -265,4 +267,36 @@ describe('subset', function () {
     const expression = math.parse('subset([1],index(0,0))')
     assert.strictEqual(expression.toTex(), '\\mathrm{subset}\\left(\\begin{bmatrix}1\\end{bmatrix},\\mathrm{index}\\left(0,0\\right)\\right)')
   })
+
+  it('should work with config legacySubset during deprecation', function () {
+    const math2 = math.create()
+    // Add a spy to temporarily disable console.warn
+    const warnStub = sinon.stub(console, 'warn')
+
+    math2.config({ legacySubset: true })
+
+    // Test legacy syntax for getting a subset of a matrix
+    const A = math2.matrix([[1, 2, 3], [4, 5, 6]])
+    const index = math2.index
+    assert.deepStrictEqual(math2.subset(A, index(1, 2)), 6)
+    assert.deepStrictEqual(math2.subset(A, index([1], 2)), 6)
+    assert.deepStrictEqual(math2.subset(A, index(1, [2])), 6)
+    assert.deepStrictEqual(math2.subset(A, index([1], [2])), 6)
+    assert.deepStrictEqual(math2.subset(A, index(1, [1, 2])).toArray(), [[5, 6]])
+    assert.deepStrictEqual(math2.subset(A, index([0, 1], 1)).toArray(), [[2], [5]])
+
+    math2.config({ legacySubset: false })
+    // Test without legacy syntax
+    assert.deepStrictEqual(math2.subset(A, index(1, 2)), 6)
+    assert.deepStrictEqual(math2.subset(A, index([1], 2)).toArray(), [6])
+    assert.deepStrictEqual(math2.subset(A, index(1, [2])).toArray(), [6])
+    assert.deepStrictEqual(math2.subset(A, index([1], [2])).toArray(), [[6]])
+    assert.deepStrictEqual(math2.subset(A, index(1, [1, 2])).toArray(), [5, 6])
+    assert.deepStrictEqual(math2.subset(A, index([1], [1, 2])).toArray(), [[5, 6]])
+    assert.deepStrictEqual(math2.subset(A, index([0, 1], 1)).toArray(), [2, 5])
+    assert.deepStrictEqual(math2.subset(A, index([0, 1], [1])).toArray(), [[2], [5]])
+
+    // Restore console.warn
+    warnStub.restore()
+  })
 })
diff --git a/test/unit-tests/json/replacer.test.js b/test/unit-tests/json/replacer.test.js
index 656b00b3f6..c918008130 100644
--- a/test/unit-tests/json/replacer.test.js
+++ b/test/unit-tests/json/replacer.test.js
@@ -59,7 +59,7 @@ describe('replacer', function () {
     const i = new math.Index(new math.Range(0, 10), 2)
     const json = '{"mathjs":"Index","dimensions":[' +
         '{"mathjs":"Range","start":0,"end":10,"step":1},' +
-        '{"mathjs":"ImmutableDenseMatrix","data":[2],"size":[1]}' +
+        '2' +
         ']}'
     assert.deepStrictEqual(JSON.stringify(i), json)
     assert.deepStrictEqual(JSON.stringify(i, replacer), json)
diff --git a/test/unit-tests/type/matrix/DenseMatrix.test.js b/test/unit-tests/type/matrix/DenseMatrix.test.js
index af2b7eded0..1540eb3146 100644
--- a/test/unit-tests/type/matrix/DenseMatrix.test.js
+++ b/test/unit-tests/type/matrix/DenseMatrix.test.js
@@ -494,10 +494,10 @@ describe('DenseMatrix', function () {
       assert.deepStrictEqual(m.size(), [3, 3])
       assert.deepStrictEqual(m.subset(index(1, 1)), 5)
       assert.deepStrictEqual(m.subset(index(new Range(0, 2), new Range(0, 2))).valueOf(), [[1, 2], [4, 5]])
-      assert.deepStrictEqual(m.subset(index(1, new Range(1, 3))).valueOf(), [[5, 6]])
-      assert.deepStrictEqual(m.subset(index(0, new Range(1, 3))).valueOf(), [[2, 3]])
-      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 1)).valueOf(), [[5], [8]])
-      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 2)).valueOf(), [[6], [9]])
+      assert.deepStrictEqual(m.subset(index(1, new Range(1, 3))).valueOf(), [5, 6])
+      assert.deepStrictEqual(m.subset(index(0, new Range(1, 3))).valueOf(), [2, 3])
+      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 1)).valueOf(), [5, 8])
+      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 2)).valueOf(), [6, 9])
       assert.deepStrictEqual(m.subset(index([0, 1, 2], [1])).valueOf(), [[2], [5], [8]])
 
       // get n-dimensional
@@ -506,9 +506,9 @@ describe('DenseMatrix', function () {
       assert.deepStrictEqual(m.subset(index(new Range(0, 2), new Range(0, 2), new Range(0, 2))).valueOf(), m.valueOf())
       assert.deepStrictEqual(m.subset(index(0, 0, 0)), 1)
       assert.deepStrictEqual(m.subset(index(1, 1, 1)).valueOf(), 8)
-      assert.deepStrictEqual(m.subset(index(1, 1, new Range(0, 2))).valueOf(), [[[7, 8]]])
-      assert.deepStrictEqual(m.subset(index(1, new Range(0, 2), 1)).valueOf(), [[[6], [8]]])
-      assert.deepStrictEqual(m.subset(index(new Range(0, 2), 1, 1)).valueOf(), [[[4]], [[8]]])
+      assert.deepStrictEqual(m.subset(index(1, 1, new Range(0, 2))).valueOf(), [7, 8])
+      assert.deepStrictEqual(m.subset(index(1, new Range(0, 2), 1)).valueOf(), [6, 8])
+      assert.deepStrictEqual(m.subset(index(new Range(0, 2), 1, 1)).valueOf(), [4, 8])
     })
 
     it('should squeeze the output when index contains a scalar', function () {
@@ -518,8 +518,8 @@ describe('DenseMatrix', function () {
 
       m = new DenseMatrix([[1, 2], [3, 4]])
       assert.deepStrictEqual(m.subset(index(1, 1)), 4)
-      assert.deepStrictEqual(m.subset(index(new Range(1, 2), 1)), new DenseMatrix([[4]]))
-      assert.deepStrictEqual(m.subset(index(1, new Range(1, 2))), new DenseMatrix([[4]]))
+      assert.deepStrictEqual(m.subset(index(new Range(1, 2), 1)), new DenseMatrix([4]))
+      assert.deepStrictEqual(m.subset(index(1, new Range(1, 2))), new DenseMatrix([4]))
       assert.deepStrictEqual(m.subset(index(new Range(1, 2), new Range(1, 2))), new DenseMatrix([[4]]))
     })
 
diff --git a/test/unit-tests/type/matrix/ImmutableDenseMatrix.test.js b/test/unit-tests/type/matrix/ImmutableDenseMatrix.test.js
index db910c47fd..57d8d33436 100644
--- a/test/unit-tests/type/matrix/ImmutableDenseMatrix.test.js
+++ b/test/unit-tests/type/matrix/ImmutableDenseMatrix.test.js
@@ -255,10 +255,10 @@ describe('ImmutableDenseMatrix', function () {
       assert.deepStrictEqual(m.size(), [3, 3])
       assert.deepStrictEqual(m.subset(index(1, 1)), 5)
       assert.deepStrictEqual(m.subset(index(new Range(0, 2), new Range(0, 2))).valueOf(), [[1, 2], [4, 5]])
-      assert.deepStrictEqual(m.subset(index(1, new Range(1, 3))).valueOf(), [[5, 6]])
-      assert.deepStrictEqual(m.subset(index(0, new Range(1, 3))).valueOf(), [[2, 3]])
-      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 1)).valueOf(), [[5], [8]])
-      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 2)).valueOf(), [[6], [9]])
+      assert.deepStrictEqual(m.subset(index(1, new Range(1, 3))).valueOf(), [5, 6])
+      assert.deepStrictEqual(m.subset(index(0, new Range(1, 3))).valueOf(), [2, 3])
+      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 1)).valueOf(), [5, 8])
+      assert.deepStrictEqual(m.subset(index(new Range(1, 3), 2)).valueOf(), [6, 9])
 
       // get n-dimensional
       m = new ImmutableDenseMatrix([[[1, 2], [3, 4]], [[5, 6], [7, 8]]])
@@ -266,9 +266,9 @@ describe('ImmutableDenseMatrix', function () {
       assert.deepStrictEqual(m.subset(index(new Range(0, 2), new Range(0, 2), new Range(0, 2))).valueOf(), m.valueOf())
       assert.deepStrictEqual(m.subset(index(0, 0, 0)), 1)
       assert.deepStrictEqual(m.subset(index(1, 1, 1)).valueOf(), 8)
-      assert.deepStrictEqual(m.subset(index(1, 1, new Range(0, 2))).valueOf(), [[[7, 8]]])
-      assert.deepStrictEqual(m.subset(index(1, new Range(0, 2), 1)).valueOf(), [[[6], [8]]])
-      assert.deepStrictEqual(m.subset(index(new Range(0, 2), 1, 1)).valueOf(), [[[4]], [[8]]])
+      assert.deepStrictEqual(m.subset(index(1, 1, new Range(0, 2))).valueOf(), [7, 8])
+      assert.deepStrictEqual(m.subset(index(1, new Range(0, 2), 1)).valueOf(), [6, 8])
+      assert.deepStrictEqual(m.subset(index(new Range(0, 2), 1, 1)).valueOf(), [4, 8])
     })
 
     it('should squeeze the output when index contains a scalar', function () {
@@ -278,8 +278,8 @@ describe('ImmutableDenseMatrix', function () {
 
       m = new ImmutableDenseMatrix([[1, 2], [3, 4]])
       assert.deepStrictEqual(m.subset(index(1, 1)), 4)
-      assert.deepStrictEqual(m.subset(index(new Range(1, 2), 1)), new ImmutableDenseMatrix([[4]]))
-      assert.deepStrictEqual(m.subset(index(1, new Range(1, 2))), new ImmutableDenseMatrix([[4]]))
+      assert.deepStrictEqual(m.subset(index(new Range(1, 2), 1)), new ImmutableDenseMatrix([4]))
+      assert.deepStrictEqual(m.subset(index(1, new Range(1, 2))), new ImmutableDenseMatrix([4]))
       assert.deepStrictEqual(m.subset(index(new Range(1, 2), new Range(1, 2))), new ImmutableDenseMatrix([[4]]))
     })
 
diff --git a/test/unit-tests/type/matrix/Index.test.js b/test/unit-tests/type/matrix/Index.test.js
index 2efef5b894..0d9ce73acd 100644
--- a/test/unit-tests/type/matrix/Index.test.js
+++ b/test/unit-tests/type/matrix/Index.test.js
@@ -8,7 +8,7 @@ const ImmutableDenseMatrix = math.ImmutableDenseMatrix
 
 describe('Index', function () {
   it('should create an Index', function () {
-    assert.deepStrictEqual(new Index(0, 2)._dimensions, [new ImmutableDenseMatrix([0]), new ImmutableDenseMatrix([2])])
+    assert.deepStrictEqual(new Index(0, 2)._dimensions, [0, 2])
 
     assert.deepStrictEqual(new Index(new Range(0, 10))._dimensions, [new Range(0, 10, 1)])
     assert.deepStrictEqual(new Index(new Range(0, 10, 2))._dimensions, [new Range(0, 10, 2)])
@@ -18,11 +18,11 @@ describe('Index', function () {
     ])
     assert.deepStrictEqual(new Index(new ImmutableDenseMatrix([0, 10]))._dimensions, [new ImmutableDenseMatrix([0, 10])])
     assert.deepStrictEqual(new Index([0, 10])._dimensions, [new ImmutableDenseMatrix([0, 10])])
-    assert.deepStrictEqual(new Index(10)._dimensions, [new ImmutableDenseMatrix([10])])
+    assert.deepStrictEqual(new Index(10)._dimensions, [10])
   })
 
   it('should create an Index from bigints', function () {
-    assert.deepStrictEqual(new Index(0n, 2n)._dimensions, [new ImmutableDenseMatrix([0]), new ImmutableDenseMatrix([2])])
+    assert.deepStrictEqual(new Index(0n, 2n)._dimensions, [0, 2])
 
     assert.deepStrictEqual(new Index(new Range(0n, 10n))._dimensions, [new Range(0, 10, 1)])
     assert.deepStrictEqual(new Index(new Range(0n, 10n, 2))._dimensions, [new Range(0, 10, 2)])
@@ -105,11 +105,11 @@ describe('Index', function () {
     assert.strictEqual(new Index(2, 5, 2).isScalar(), true)
     assert.strictEqual(new Index(2).isScalar(), true)
     assert.strictEqual(new Index([0, 1, 2], 2).isScalar(), false)
-    assert.strictEqual(new Index([3], [2]).isScalar(), true)
+    assert.strictEqual(new Index([3], [2]).isScalar(), false)
     assert.strictEqual(new Index([0, 1, 2], [2]).isScalar(), false)
     assert.strictEqual(new Index(new Range(2, 10)).isScalar(), false)
     assert.strictEqual(new Index(new ImmutableDenseMatrix([2, 10])).isScalar(), false)
-    assert.strictEqual(new Index(new ImmutableDenseMatrix([2])).isScalar(), true)
+    assert.strictEqual(new Index(new ImmutableDenseMatrix([2])).isScalar(), false)
     assert.strictEqual(new Index(2, new Range(0, 4), 2).isScalar(), false)
     assert.strictEqual(new Index(2, new ImmutableDenseMatrix([0, 4]), 2).isScalar(), false)
     assert.strictEqual(new Index(new Range(0, 2), new Range(0, 4)).isScalar(), false)
@@ -119,8 +119,8 @@ describe('Index', function () {
   })
 
   it('should clone an Index', function () {
-    const index1 = new Index(2, new Range(0, 4), new ImmutableDenseMatrix([0, 2]))
-    const index2 = index1.clone(0)
+    const index1 = new Index([2], new Range(0, 4), new ImmutableDenseMatrix([0, 2]))
+    const index2 = index1.clone()
 
     assert.deepStrictEqual(index1, index2)
     assert.notStrictEqual(index1, index2)
@@ -131,9 +131,9 @@ describe('Index', function () {
 
   it('should stringify an index', function () {
     assert.strictEqual(new Index().toString(), '[]')
-    assert.strictEqual(new Index(2, 3).toString(), '[[2], [3]]')
-    assert.strictEqual(new Index(2, 3, 1).toString(), '[[2], [3], [1]]')
-    assert.strictEqual(new Index(2, new Range(0, 3)).toString(), '[[2], 0:3]')
+    assert.strictEqual(new Index(2, 3).toString(), '[2, 3]')
+    assert.strictEqual(new Index(2, 3, 1).toString(), '[2, 3, 1]')
+    assert.strictEqual(new Index(2, new Range(0, 3)).toString(), '[2, 0:3]')
     assert.strictEqual(new Index(new Range(0, 6, 2)).toString(), '[0:2:6]')
     assert.strictEqual(new Index(new ImmutableDenseMatrix([0, 6, 2])).toString(), '[[0, 6, 2]]')
     assert.deepStrictEqual(new Index('property').toString(), '["property"]')
@@ -145,7 +145,7 @@ describe('Index', function () {
         mathjs: 'Index',
         dimensions: [
           new Range(0, 10, 1),
-          new ImmutableDenseMatrix([2]),
+          2,
           new ImmutableDenseMatrix([1, 2, 3])
         ]
       })
@@ -155,7 +155,7 @@ describe('Index', function () {
     const json = {
       dimensions: [
         new Range(0, 10, 1),
-        new ImmutableDenseMatrix([2]),
+        2,
         new ImmutableDenseMatrix([1, 2, 3])
       ]
     }
@@ -169,8 +169,8 @@ describe('Index', function () {
   it('should get the range for a given dimension', function () {
     const index = new Index(2, new Range(0, 8, 2), new Range(3, -1, -1), new ImmutableDenseMatrix([1, 2, 3]))
 
-    assert(index.dimension(0) instanceof ImmutableDenseMatrix)
-    assert.deepStrictEqual(index.dimension(0), new ImmutableDenseMatrix([2]))
+    assert(Number.isInteger(index.dimension(0)))
+    assert.deepStrictEqual(index.dimension(0), 2)
 
     assert(index.dimension(1) instanceof Range)
     assert.deepStrictEqual(index.dimension(1), new Range(0, 8, 2))
@@ -197,7 +197,7 @@ describe('Index', function () {
     })
 
     assert.deepStrictEqual(log, [
-      { dimension: new ImmutableDenseMatrix([2]), index: 0 },
+      { dimension: 2, index: 0 },
       { dimension: new Range(0, 8, 2), index: 1 },
       { dimension: new Range(3, -1, -1), index: 2 },
       { dimension: new ImmutableDenseMatrix([1, 2, 3]), index: 3 }
@@ -229,12 +229,12 @@ describe('Index', function () {
     assert.deepStrictEqual(new Index(new Range(2, 5), new Range(0, 8, 2), 2, new ImmutableDenseMatrix([1, 2])).toArray(), [
       [2, 3, 4],
       [0, 2, 4, 6],
-      [2],
+      2,
       [1, 2]
     ])
 
     assert.deepStrictEqual(new Index(2, new Range(0, 8, 2), new Range(3, -1, -1), new Range(2, 4, 0), [1, 2]).toArray(), [
-      [2],
+      2,
       [0, 2, 4, 6],
       [3, 2, 1, 0],
       [],
@@ -246,7 +246,7 @@ describe('Index', function () {
 
   it('valueOf should return the expanded array', function () {
     assert.deepStrictEqual(new Index(2, new Range(0, 8, 2), new Range(3, -1, -1), [1, 2], new ImmutableDenseMatrix([3, 4])).valueOf(), [
-      [2],
+      2,
       [0, 2, 4, 6],
       [3, 2, 1, 0],
       [1, 2],
diff --git a/test/unit-tests/type/matrix/function/index.test.js b/test/unit-tests/type/matrix/function/index.test.js
index 3ef9c4221e..b2b19bb726 100644
--- a/test/unit-tests/type/matrix/function/index.test.js
+++ b/test/unit-tests/type/matrix/function/index.test.js
@@ -3,7 +3,6 @@ import assert from 'assert'
 
 import math from '../../../../../src/defaultInstance.js'
 const Range = math.Range
-const ImmutableDenseMatrix = math.ImmutableDenseMatrix
 
 describe('index', function () {
   it('should create an index', function () {
@@ -22,13 +21,13 @@ describe('index', function () {
   it('should create an index from bignumbers (downgrades to numbers)', function () {
     const index = math.index(new Range(math.bignumber(2), math.bignumber(6)), math.bignumber(3))
     assert.ok(index instanceof math.Index)
-    assert.deepStrictEqual(index._dimensions, [new Range(2, 6, 1), new ImmutableDenseMatrix([3])])
+    assert.deepStrictEqual(index._dimensions, [new Range(2, 6, 1), 3])
   })
 
   it('should create an index from bigints (downgrades to numbers)', function () {
     const index = math.index(new Range(2n, 6n), 3n)
     assert.ok(index instanceof math.Index)
-    assert.deepStrictEqual(index._dimensions, [new Range(2, 6, 1), new ImmutableDenseMatrix([3])])
+    assert.deepStrictEqual(index._dimensions, [new Range(2, 6, 1), 3])
   })
 
   it('should LaTeX index', function () {