diff --git a/arrays.scad b/arrays.scad
index d5e03a5..6e57953 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -1131,7 +1131,8 @@ function add_scalar(v,s) =
 // Description:
 //   Extracts the entries listed in idx from each entry in M.  For a matrix this means
 //   selecting a specified set of columsn.  If idx is a number the return is a vector, otherwise
-//   it is a list of lists (the submatrix).  
+//   it is a list of lists (the submatrix).  Note that unlike subindex, even if you give a number for
+//   an index the output includes all levels of list nesting.  
 // Arguments:
 //   M = The given list of lists.
 //   idx = The index, list of indices, or range of indices to fetch.
@@ -1146,6 +1147,36 @@ function subindex(M, idx) =
       ? [for(row=M) row[idx]]
       : [for(row=M) [for(i=idx) row[i]]];
 
+
+// Function: submatrix()
+// Usage: submatrix(M, idx1, idx2)
+// Description:
+//   The input must be a list of lists (a matrix or 2d array).  Returns a submatrix by selecting the rows listed in idx1 and columsn listed in idx2.
+// Arguments:
+//   M = Given list of lists
+//   idx1 = rows index list or range
+//   idx2 = column index list or range
+// Example:
+//   M = [[ 1, 2, 3, 4, 5],
+//        [ 6, 7, 8, 9,10],
+//        [11,12,13,14,15],
+//        [16,17,18,19,20],
+//        [21,22,23,24,25]];
+//   submatrix(M,[1:2],[3:4]);  // Returns [[9, 10], [14, 15]]
+//   submatrix(M,[1], [3,4]));  // Returns [[9,10]]
+//   submatrix(M,1, [3,4]));  // Returns [[9,10]]
+//   submatrix(M,1,3));  // Returns [[9]]
+//   submatrix(M, [3,4],1); // Returns  [[17],[22]]);
+//   submatrix(M, [1,3],[2,4]); // Returns [[8,10],[18,20]]);
+//   A = [[true,    17, "test"],
+//        [[4,2],   91, false],
+//        [6,    [3,4], undef]];
+//   submatrix(A,[0,2],[1,2]);   // Returns [[17, "test"], [[3, 4], undef]]
+
+function submatrix(M,idx1,idx2) =
+    [for(i=idx1) [for(j=idx2) M[i][j] ] ];
+
+
 // Function: zip()
 // Usage:
 //   zip(v1, v2, v3, [fit], [fill]);
diff --git a/math.scad b/math.scad
index 7e16bbe..e40c1c6 100644
--- a/math.scad
+++ b/math.scad
@@ -716,18 +716,6 @@ function matrix_inverse(A) =
     linear_solve(A,ident(len(A)));
 
 
-// Function: submatrix()
-// Usage: submatrix(M, ind1, ind2)
-// Description:
-//   Returns a submatrix with the specified index ranges or index sets.  
-function submatrix(M,ind1,ind2) =
-    assert( is_matrix(M), "Input must be a matrix." )
-    [for(i=ind1) 
-        [for(j=ind2) 
-            assert( ! is_undef(M[i][j]), "Invalid indexing." )
-            M[i][j] ] ];
-    
-
 // Function: qr_factor()
 // Usage: qr = qr_factor(A)
 // Description:
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 32bc5ff..64a7fbb 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -357,6 +357,27 @@ module test_subindex() {
 test_subindex();
 
 
+// Need decision about behavior for out of bounds ranges, empty ranges
+module test_submatrix(){
+  M = [[1,2,3,4,5],
+       [6,7,8,9,10],
+       [11,12,13,14,15],
+       [16,17,18,19,20],
+       [21,22,23,24,25]];
+  assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]);
+  assert_equal(submatrix(M,[1], [3,4]), [[9,10]]);
+  assert_equal(submatrix(M,1, [3,4]), [[9,10]]);
+  assert_equal(submatrix(M, [3,4],1), [[17],[22]]);
+  assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]);
+  assert_equal(submatrix(M, 1,3), [[9]]);
+  A = [[true,    17, "test"],
+     [[4,2],   91, false],
+     [6,    [3,4], undef]];
+  assert_equal(submatrix(A,[0,2],[1,2]),[[17, "test"], [[3, 4], undef]]);
+}
+test_submatrix();
+
+
 module test_force_list() {
     assert_equal(force_list([3,4,5]), [3,4,5]);
     assert_equal(force_list(5), [5]);
diff --git a/tests/test_math.scad b/tests/test_math.scad
index 91e662f..12f7d2e 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -853,22 +853,6 @@ module test_real_roots(){
 }
 test_real_roots();
 
-// Need decision about behavior for out of bounds ranges, empty ranges
-module test_submatrix(){
-  M = [[1,2,3,4,5],
-       [6,7,8,9,10],
-       [11,12,13,14,15],
-       [16,17,18,19,20],
-       [21,22,23,24,25]];
-  assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]);
-  assert_equal(submatrix(M,[1], [3,4]), [[9,10]]);
-  assert_equal(submatrix(M,1, [3,4]), [[9,10]]);
-  assert_equal(submatrix(M, [3,4],1), [[17],[22]]);
-  assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]);
-}
-test_submatrix();
-
-
 
 module test_qr_factor() {
   // Check that R is upper triangular