diff --git a/src/fastmath/matrix.clj b/src/fastmath/matrix.clj
index df517e7d..e3497e54 100644
--- a/src/fastmath/matrix.clj
+++ b/src/fastmath/matrix.clj
@@ -1217,12 +1217,14 @@
 (defn eigenvalues-matrix
   "Returns eigenvalues for given matrix as a diagonal or block diagonal matrix"
   [A]
-  (->> (mat->RealMatrix A)
-       (EigenDecomposition.)
-       ^RealMatrix (.getD)
-       (.getData)
-       (m/double-double-array->seq)
-       (apply rows->mat)))
+  (let [data (->> (mat->RealMatrix A)
+                  (EigenDecomposition.)
+                  ^RealMatrix (.getD)
+                  (.getData))
+        n (alength ^"[[D" data)]
+    (if (<= n 4)
+      (apply rows->mat (m/double-double-array->seq data))
+       data)))
 
 (defn square?
   "Is matrix square?"
@@ -1345,12 +1347,14 @@
   "Returns eigenvectors as a matrix (columns). Vectors can be normalized."
   ([A] (eigenvectors A false))
   ([A normalize?]
-   (let [evs (->> (mat->RealMatrix A)
+   (let [data (->> (mat->RealMatrix A)
                   (EigenDecomposition.)
                   ^RealMatrix (.getV)
-                  (.getData)
-                  (m/double-double-array->seq)
-                  (apply rows->mat))]
+                  (.getData))
+          n (alength ^"[[D" data)
+          evs (if (<= n 4)
+                  (apply rows->mat (m/double-double-array->seq data))
+                  data)]
      (if normalize? (normalize evs) evs))))
 
 ;;
diff --git a/test/fastmath/matrix_test.clj b/test/fastmath/matrix_test.clj
index 9738d29d..69d7203f 100644
--- a/test/fastmath/matrix_test.clj
+++ b/test/fastmath/matrix_test.clj
@@ -360,3 +360,33 @@
     :frobenius m44ra 12.12436
     :max m44ra 5
     [1] m44ra 20.32997715038671))
+
+;;
+
+(t/deftest eigenvalues-and-eigenvectors
+  ;; Test small matrices (2x2, 3x3, 4x4) - existing behavior
+  (t/are [m] (let [evecs (sut/eigenvectors m)
+                   evals (sut/eigenvalues m)
+                   D (sut/eigenvalues-matrix m)]
+               (and (= (sut/nrow m) (sut/nrow evecs) (sut/ncol evecs))
+                    (= (sut/nrow m) (sut/nrow D) (sut/ncol D))
+                    (= (count evals) (sut/nrow m))))
+    m22 m33 m44 m44a m44ra)
+
+  ;; Test matrices larger than 4x4 (issue #42)
+  (let [make-diagonal (fn [n]
+                        (let [arr (make-array Double/TYPE n n)]
+                          (doseq [i (range n)]
+                            (aset arr i i (double (inc i))))
+                          arr))]
+    (doseq [n [5 6 10]]
+      (let [diag-mat (make-diagonal n)
+            evecs (sut/eigenvectors diag-mat)
+            evals (sut/eigenvalues diag-mat)
+            D (sut/eigenvalues-matrix diag-mat)]
+        ;; Should not throw, dimensions should be correct
+        (t/is (= n (sut/nrow evecs)) (str "eigenvectors nrow for " n "x" n))
+        (t/is (= n (sut/ncol evecs)) (str "eigenvectors ncol for " n "x" n))
+        (t/is (= n (sut/nrow D)) (str "eigenvalues-matrix nrow for " n "x" n))
+        (t/is (= n (sut/ncol D)) (str "eigenvalues-matrix ncol for " n "x" n))
+        (t/is (= n (count evals)) (str "eigenvalues count for " n "x" n))))))