diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml
index 67f782048..3b04949b9 100644
--- a/.github/workflows/python-package-conda.yml
+++ b/.github/workflows/python-package-conda.yml
@@ -49,6 +49,9 @@ jobs:
         # Use xvfb-run instead of pytest-xvfb to get proper mpl backend
         # Use coverage instead of pytest-cov to get .coverage file
         # See https://github.com/python-control/python-control/pull/504
+       
+        python3 -c 'import numpy; numpy.show_config()'
+
         xvfb-run --auto-servernum coverage run -m pytest control/tests
 
     - name: Coveralls parallel
diff --git a/control/statesp.py b/control/statesp.py
index f04bc55b3..09abbb706 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -887,6 +887,8 @@ def horner(self, x, warn_infinite=True):
         """
         # Make sure the argument is a 1D array of complex numbers
         x_arr = np.atleast_1d(x).astype(complex, copy=False)
+        if len(x_arr.shape) > 1:
+            raise ValueError("input list must be 1D")
 
         # return fast on systems with 0 or 1 state
         if not config.defaults['statesp.use_numpy_matrix']:
@@ -908,10 +910,6 @@ def horner(self, x, warn_infinite=True):
             # Fall back because either Slycot unavailable or cannot handle
             # certain cases.
 
-            # Make sure that we are operating on a simple list
-            if len(x_arr.shape) > 1:
-                raise ValueError("input list must be 1D")
-
             # Preallocate
             out = empty((self.noutputs, self.ninputs, len(x_arr)),
                         dtype=complex)
@@ -919,10 +917,14 @@ def horner(self, x, warn_infinite=True):
             # TODO: can this be vectorized?
             for idx, x_idx in enumerate(x_arr):
                 try:
-                    out[:, :, idx] = np.dot(
-                        self.C,
-                        solve(x_idx * eye(self.nstates) - self.A, self.B)) \
-                        + self.D
+                    xI_A = x_idx * eye(self.nstates) - self.A
+                    xI_A_inv = solve(xI_A, self.B)
+                    # gh-664: xI_A did not raise singular matrix error,
+                    print(f"DEBUG np.linalg.solve(\n,"
+                          f"{xI_A}\n"
+                          f"{self.B} = \n"
+                          f"{xI_A_inv}")
+                    out[:, :, idx] = np.dot(self.C, xI_A_inv) + self.D
                 except LinAlgError:
                     # Issue a warning messsage, for consistency with xferfcn
                     if warn_infinite: