diff --git a/control/lti.py b/control/lti.py
index 445775f5f..30569863a 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -644,8 +644,10 @@ def dcgain(sys):
     Returns
     -------
     gain : ndarray
-        The zero-frequency gain, or np.nan if the system has a pole
-        at the origin
+        The zero-frequency gain, or (inf + nanj) if the system has a pole at
+        the origin, (nan + nanj) if there is a pole/zero cancellation at the
+        origin.
+
     """
     return sys.dcgain()
 
diff --git a/control/margins.py b/control/margins.py
index 02d615e05..1231c5388 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -294,8 +294,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
             num_iw, den_iw = _poly_iw(sys)
             # frequency for gain margin: phase crosses -180 degrees
             w_180 = _poly_iw_real_crossing(num_iw, den_iw, epsw)
-            with np.errstate(all='ignore'):  # den=0 is okay
-                w180_resp = sys(1J * w_180)
+            w180_resp = sys(1J * w_180, warn_infinite=False)  # den=0 is okay
 
             # frequency for phase margin : gain crosses magnitude 1
             wc = _poly_iw_mag1_crossing(num_iw, den_iw, epsw)
diff --git a/control/statesp.py b/control/statesp.py
index abd55ad15..d2b613024 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -668,7 +668,7 @@ def __rdiv__(self, other):
         raise NotImplementedError(
             "StateSpace.__rdiv__ is not implemented yet.")
 
-    def __call__(self, x, squeeze=None):
+    def __call__(self, x, squeeze=None, warn_infinite=True):
         """Evaluate system's transfer function at complex frequency.
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
@@ -689,6 +689,8 @@ def __call__(self, x, squeeze=None):
             keep all indices (output, input and, if omega is array_like,
             frequency) even if the system is SISO. The default value can be
             set using config.defaults['control.squeeze_frequency_response'].
+        warn_infinite : bool, optional
+            If set to `False`, don't warn if frequency response is infinite.
 
         Returns
         -------
@@ -702,7 +704,7 @@ def __call__(self, x, squeeze=None):
 
         """
         # Use Slycot if available
-        out = self.horner(x)
+        out = self.horner(x, warn_infinite=warn_infinite)
         return _process_frequency_response(self, x, out, squeeze=squeeze)
 
     def slycot_laub(self, x):
@@ -723,7 +725,9 @@ def slycot_laub(self, x):
             Frequency response
         """
         from slycot import tb05ad
-        x_arr = np.atleast_1d(x) # array-like version of x
+
+        # Make sure the argument is a 1D array of complex numbers
+        x_arr = np.atleast_1d(x).astype(complex, copy=False)
 
         # Make sure that we are operating on a simple list
         if len(x_arr.shape) > 1:
@@ -758,7 +762,7 @@ def slycot_laub(self, x):
             out[:, :, kk+1] = result[0] + self.D
         return out
 
-    def horner(self, x):
+    def horner(self, x, warn_infinite=True):
         """Evaluate system's transfer function at complex frequency
         using Laub's or Horner's method.
 
@@ -788,21 +792,38 @@ def horner(self, x):
         except (ImportError, Exception):
             # Fall back because either Slycot unavailable or cannot handle
             # certain cases.
-            x_arr = np.atleast_1d(x) # force to be an array
+
+            # Make sure the argument is a 1D array of complex numbers
+            x_arr = np.atleast_1d(x).astype(complex, copy=False)
 
             # 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)
+            out = empty((self.noutputs, self.ninputs, len(x_arr)),
+                        dtype=complex)
 
             #TODO: can this be vectorized?
             for idx, x_idx in enumerate(x_arr):
-                out[:,:,idx] = \
-                    np.dot(self.C,
+                try:
+                    out[:,:,idx] = np.dot(
+                        self.C,
                         solve(x_idx * eye(self.nstates) - self.A, self.B)) \
-                    + self.D
+                        + self.D
+                except LinAlgError:
+                    # Issue a warning messsage, for consistency with xferfcn
+                    if warn_infinite:
+                        warn("singular matrix in frequency response",
+                             RuntimeWarning)
+
+                    # Evaluating at a pole.  Return value depends if there
+                    # is a zero at the same point or not.
+                    if x_idx in self.zero():
+                        out[:,:,idx] = complex(np.nan, np.nan)
+                    else:
+                        out[:,:,idx] = complex(np.inf, np.nan)
+
         return out
 
     def freqresp(self, omega):
@@ -1183,7 +1204,7 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
         Ad, Bd, C, D, _ = cont2discrete(sys, Twarp, method, alpha)
         return StateSpace(Ad, Bd, C, D, Ts)
 
-    def dcgain(self):
+    def dcgain(self, warn_infinite=False):
         """Return the zero-frequency gain
 
         The zero-frequency gain of a continuous-time state-space
@@ -1198,20 +1219,13 @@ def dcgain(self):
         Returns
         -------
         gain : ndarray
-            An array of shape (outputs,inputs); the array will either
-            be the zero-frequency (or DC) gain, or, if the frequency
-            response is singular, the array will be filled with np.nan.
+            An array of shape (outputs,inputs); the array will either be the
+            zero-frequency (or DC) gain, or, if the frequency response is
+            singular, the array will be filled with (inf + nanj).
+
         """
-        try:
-            if self.isctime():
-                gain = np.asarray(self.D -
-                                  self.C.dot(np.linalg.solve(self.A, self.B)))
-            else:
-                gain = np.squeeze(self.horner(1))
-        except LinAlgError:
-            # eigenvalue at DC
-            gain = np.tile(np.nan, (self.noutputs, self.ninputs))
-        return np.squeeze(gain)
+        return self(0, warn_infinite=warn_infinite) if self.isctime() \
+            else self(1, warn_infinite=warn_infinite)
 
     def _isstatic(self):
         """True if and only if the system has no dynamics, that is,
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index 86de0e77a..54fe3b0e0 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -364,3 +364,140 @@ def test_phase_wrap(TF, wrap_phase, min_phase, max_phase):
     mag, phase, omega = ctrl.bode(TF, wrap_phase=wrap_phase)
     assert(min(phase) >= min_phase)
     assert(max(phase) <= max_phase)
+
+
+def test_freqresp_warn_infinite():
+    """Test evaluation warnings for transfer functions w/ pole at the origin"""
+    sys_finite = ctrl.tf([1], [1, 0.01])
+    sys_infinite = ctrl.tf([1], [1, 0.01, 0])
+
+    # Transfer function with finite zero frequency gain
+    np.testing.assert_almost_equal(sys_finite(0), 100)
+    np.testing.assert_almost_equal(sys_finite(0, warn_infinite=False), 100)
+    np.testing.assert_almost_equal(sys_finite(0, warn_infinite=True), 100)
+
+    # Transfer function with infinite zero frequency gain
+    with pytest.warns(RuntimeWarning, match="divide by zero"):
+        np.testing.assert_almost_equal(
+            sys_infinite(0), complex(np.inf, np.nan))
+    with pytest.warns(RuntimeWarning, match="divide by zero"):
+        np.testing.assert_almost_equal(
+            sys_infinite(0, warn_infinite=True), complex(np.inf, np.nan))
+    np.testing.assert_almost_equal(
+        sys_infinite(0, warn_infinite=False), complex(np.inf, np.nan))
+
+    # Switch to state space
+    sys_finite = ctrl.tf2ss(sys_finite)
+    sys_infinite = ctrl.tf2ss(sys_infinite)
+
+    # State space system with finite zero frequency gain
+    np.testing.assert_almost_equal(sys_finite(0), 100)
+    np.testing.assert_almost_equal(sys_finite(0, warn_infinite=False), 100)
+    np.testing.assert_almost_equal(sys_finite(0, warn_infinite=True), 100)
+
+    # State space system with infinite zero frequency gain
+    with pytest.warns(RuntimeWarning, match="singular matrix"):
+        np.testing.assert_almost_equal(
+            sys_infinite(0), complex(np.inf, np.nan))
+    with pytest.warns(RuntimeWarning, match="singular matrix"):
+        np.testing.assert_almost_equal(
+            sys_infinite(0, warn_infinite=True), complex(np.inf, np.nan))
+    np.testing.assert_almost_equal(sys_infinite(
+        0, warn_infinite=False), complex(np.inf, np.nan))
+
+
+def test_dcgain_consistency():
+    """Test to make sure that DC gain is consistently evaluated"""
+    # Set up transfer function with pole at the origin
+    sys_tf = ctrl.tf([1], [1, 0])
+    assert 0 in sys_tf.pole()
+
+    # Set up state space system with pole at the origin
+    sys_ss = ctrl.tf2ss(sys_tf)
+    assert 0 in sys_ss.pole()
+
+    # Finite (real) numerator over 0 denominator => inf + nanj
+    np.testing.assert_equal(
+        sys_tf(0, warn_infinite=False), complex(np.inf, np.nan))
+    np.testing.assert_equal(
+        sys_ss(0, warn_infinite=False), complex(np.inf, np.nan))
+    np.testing.assert_equal(
+        sys_tf(0j, warn_infinite=False), complex(np.inf, np.nan))
+    np.testing.assert_equal(
+        sys_ss(0j, warn_infinite=False), complex(np.inf, np.nan))
+    np.testing.assert_equal(
+        sys_tf.dcgain(warn_infinite=False), complex(np.inf, np.nan))
+    np.testing.assert_equal(
+        sys_ss.dcgain(warn_infinite=False), complex(np.inf, np.nan))
+
+    # Set up transfer function with pole, zero at the origin
+    sys_tf = ctrl.tf([1, 0], [1, 0])
+    assert 0 in sys_tf.pole()
+    assert 0 in sys_tf.zero()
+
+    # Pole and zero at the origin should give nan + nanj for the response
+    np.testing.assert_equal(
+        sys_tf(0, warn_infinite=False), complex(np.nan, np.nan))
+    np.testing.assert_equal(
+        sys_tf(0j, warn_infinite=False), complex(np.nan, np.nan))
+    np.testing.assert_equal(
+        sys_tf.dcgain(warn_infinite=False), complex(np.nan, np.nan))
+
+    # Set up state space version
+    sys_ss = ctrl.tf2ss(ctrl.tf([1, 0], [1, 1])) * \
+        ctrl.tf2ss(ctrl.tf([1], [1, 0]))
+
+    # Different systems give different representations => test accordingly
+    if 0 in sys_ss.pole() and 0 in sys_ss.zero():
+        # Pole and zero at the origin => should get (nan + nanj)
+        np.testing.assert_equal(
+            sys_ss(0, warn_infinite=False), complex(np.nan, np.nan))
+        np.testing.assert_equal(
+            sys_ss(0j, warn_infinite=False), complex(np.nan, np.nan))
+        np.testing.assert_equal(
+            sys_ss.dcgain(warn_infinite=False), complex(np.nan, np.nan))
+    elif 0 in sys_ss.pole():
+        # Pole at the origin, but zero elsewhere => should get (inf + nanj)
+        np.testing.assert_equal(
+            sys_ss(0, warn_infinite=False), complex(np.inf, np.nan))
+        np.testing.assert_equal(
+            sys_ss(0j, warn_infinite=False), complex(np.inf, np.nan))
+        np.testing.assert_equal(
+            sys_ss.dcgain(warn_infinite=False), complex(np.inf, np.nan))
+    else:
+        # Near pole/zero cancellation => nothing sensible to check
+        pass
+
+    # Pole with non-zero, complex numerator => inf + infj
+    s = ctrl.tf('s')
+    sys_tf = (s + 1) / (s**2 + 1)
+    assert 1j in sys_tf.pole()
+
+    # Set up state space system with pole on imaginary axis
+    sys_ss = ctrl.tf2ss(sys_tf)
+    assert 1j in sys_tf.pole()
+
+    # Make sure we get correct response if evaluated at the pole
+    np.testing.assert_equal(
+        sys_tf(1j, warn_infinite=False), complex(np.inf, np.inf))
+
+    # For state space, numerical errors come into play
+    resp_ss = sys_ss(1j, warn_infinite=False)
+    if np.isfinite(resp_ss):
+        assert abs(resp_ss) > 1e15
+    else:
+        if resp_ss != complex(np.inf, np.inf):
+            pytest.xfail("statesp evaluation at poles not fully implemented")
+        else:
+            np.testing.assert_equal(resp_ss, complex(np.inf, np.inf))
+
+    # DC gain is finite
+    np.testing.assert_almost_equal(sys_tf.dcgain(), 1.)
+    np.testing.assert_almost_equal(sys_ss.dcgain(), 1.)
+
+    # Make sure that we get the *signed* DC gain
+    sys_tf = -1 / (s + 1)
+    np.testing.assert_almost_equal(sys_tf.dcgain(), -1)
+
+    sys_ss = ctrl.tf2ss(sys_tf)
+    np.testing.assert_almost_equal(sys_ss.dcgain(), -1)
diff --git a/control/tests/input_element_int_test.py b/control/tests/input_element_int_test.py
index ecfaab834..94e5efcb5 100644
--- a/control/tests/input_element_int_test.py
+++ b/control/tests/input_element_int_test.py
@@ -23,7 +23,7 @@ def test_tf_den_with_numpy_int_element(self):
 
         sys = tf(num, den)
 
-        np.testing.assert_array_max_ulp(1., dcgain(sys))
+        np.testing.assert_almost_equal(1., dcgain(sys))
 
     def test_tf_num_with_numpy_int_element(self):
         num = np.convolve([1], [1, 1])
@@ -31,7 +31,7 @@ def test_tf_num_with_numpy_int_element(self):
 
         sys = tf(num, den)
 
-        np.testing.assert_array_max_ulp(1., dcgain(sys))
+        np.testing.assert_almost_equal(1., dcgain(sys))
 
     # currently these pass
     def test_tf_input_with_int_element(self):
@@ -40,7 +40,7 @@ def test_tf_input_with_int_element(self):
 
         sys = tf(num, den)
 
-        np.testing.assert_array_max_ulp(1., dcgain(sys))
+        np.testing.assert_almost_equal(1., dcgain(sys))
 
     def test_ss_input_with_int_element(self):
         a = np.array([[0, 1],
@@ -52,7 +52,7 @@ def test_ss_input_with_int_element(self):
 
         sys = ss(a, b, c, d)
         sys2 = tf(sys)
-        np.testing.assert_array_max_ulp(dcgain(sys), dcgain(sys2))
+        np.testing.assert_almost_equal(dcgain(sys), dcgain(sys2))
 
     def test_ss_input_with_0int_dcgain(self):
         a = np.array([[0, 1],
@@ -63,4 +63,4 @@ def test_ss_input_with_0int_dcgain(self):
         d = 0
         sys = ss(a, b, c, d)
         np.testing.assert_allclose(dcgain(sys), 0,
-                                   atol=np.finfo(np.float).epsneg)
\ No newline at end of file
+                                   atol=np.finfo(np.float).epsneg)
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 1c76efbc0..983b9d7a6 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -498,7 +498,7 @@ def test_dc_gain_cont(self):
         np.testing.assert_allclose(sys2.dcgain(), expected)
 
         sys3 = StateSpace(0., 1., 1., 0.)
-        np.testing.assert_equal(sys3.dcgain(), np.nan)
+        np.testing.assert_equal(sys3.dcgain(), complex(np.inf, np.nan))
 
     def test_dc_gain_discr(self):
         """Test DC gain for discrete-time state-space systems."""
@@ -516,14 +516,14 @@ def test_dc_gain_discr(self):
 
         # summer
         sys = StateSpace(1, 1, 1, 0, True)
-        np.testing.assert_equal(sys.dcgain(), np.nan)
+        np.testing.assert_equal(sys.dcgain(), complex(np.inf, np.nan))
 
     @pytest.mark.parametrize("outputs", range(1, 6))
     @pytest.mark.parametrize("inputs", range(1, 6))
     @pytest.mark.parametrize("dt", [None, 0, 1, True],
                              ids=["dtNone", "c", "dt1", "dtTrue"])
     def test_dc_gain_integrator(self, outputs, inputs, dt):
-        """DC gain when eigenvalue at DC returns appropriately sized array of nan.
+        """DC gain w/ pole at origin returns appropriately sized array of inf.
 
         the SISO case is also tested in test_dc_gain_{cont,discr}
         time systems (dt=0)
@@ -539,8 +539,15 @@ def test_dc_gain_integrator(self, outputs, inputs, dt):
         c = np.eye(max(outputs, states))[:outputs, :states]
         d = np.zeros((outputs, inputs))
         sys = StateSpace(a, b, c, d, dt)
-        dc = np.squeeze(np.full_like(d, np.nan))
-        np.testing.assert_array_equal(dc, sys.dcgain())
+        dc = np.full_like(d, complex(np.inf, np.nan), dtype=complex)
+        if sys.issiso():
+            dc = dc.squeeze()
+
+        try:
+            np.testing.assert_array_equal(dc, sys.dcgain())
+        except NotImplementedError:
+            # Skip MIMO tests if there is no slycot
+            pytest.skip("slycot required for MIMO dcgain")
 
     def test_scalar_static_gain(self):
         """Regression: can we create a scalar static gain?
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index 782fcaa13..b892655e9 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -807,7 +807,7 @@ def test_dcgain_cont(self):
         np.testing.assert_equal(sys2.dcgain(), 2)
 
         sys3 = TransferFunction(6, [1, 0])
-        np.testing.assert_equal(sys3.dcgain(), np.inf)
+        np.testing.assert_equal(sys3.dcgain(), complex(np.inf, np.nan))
 
         num = [[[15], [21], [33]], [[10], [14], [22]]]
         den = [[[1, 3], [2, 3], [3, 3]], [[1, 5], [2, 7], [3, 11]]]
@@ -827,8 +827,13 @@ def test_dcgain_discr(self):
 
         # differencer
         sys = TransferFunction(1, [1, -1], True)
+        np.testing.assert_equal(sys.dcgain(), complex(np.inf, np.nan))
+
+        # differencer, with warning
+        sys = TransferFunction(1, [1, -1], True)
         with pytest.warns(RuntimeWarning, match="divide by zero"):
-            np.testing.assert_equal(sys.dcgain(), np.inf)
+            np.testing.assert_equal(
+                sys.dcgain(warn_infinite=True), complex(np.inf, np.nan))
 
         # summer
         sys = TransferFunction([1, -1], [1], True)
diff --git a/control/xferfcn.py b/control/xferfcn.py
index c9fbe9006..50e4870a8 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -234,7 +234,7 @@ def __init__(self, *args, **kwargs):
                     dt = config.defaults['control.default_dt']
         self.dt = dt
 
-    def __call__(self, x, squeeze=None):
+    def __call__(self, x, squeeze=None, warn_infinite=True):
         """Evaluate system's transfer function at complex frequencies.
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
@@ -262,6 +262,8 @@ def __call__(self, x, squeeze=None):
             If True and the system is single-input single-output (SISO),
             return a 1D array rather than a 3D array.  Default value (True)
             set by config.defaults['control.squeeze_frequency_response'].
+        warn_infinite : bool, optional
+            If set to `False`, turn off divide by zero warning.
 
         Returns
         -------
@@ -274,10 +276,10 @@ def __call__(self, x, squeeze=None):
             then single-dimensional axes are removed.
 
         """
-        out = self.horner(x)
+        out = self.horner(x, warn_infinite=warn_infinite)
         return _process_frequency_response(self, x, out, squeeze=squeeze)
 
-    def horner(self, x):
+    def horner(self, x, warn_infinite=True):
         """Evaluate system's transfer function at complex frequency
         using Horner's method.
 
@@ -298,17 +300,22 @@ def horner(self, x):
             Frequency response
 
         """
-        x_arr = np.atleast_1d(x)        # force to be an array
+        # Make sure the argument is a 1D array of complex numbers
+        x_arr = np.atleast_1d(x).astype(complex, copy=False)
 
         # Make sure that we are operating on a simple list
         if len(x_arr.shape) > 1:
             raise ValueError("input list must be 1D")
 
+        # Initialize the output matrix in the proper shape
         out = empty((self.noutputs, self.ninputs, len(x_arr)), dtype=complex)
-        for i in range(self.noutputs):
-            for j in range(self.ninputs):
-                out[i][j] = (polyval(self.num[i][j], x) /
-                             polyval(self.den[i][j], x))
+
+        # Set up error processing based on warn_infinite flag
+        with np.errstate(all='warn' if warn_infinite else 'ignore'):
+            for i in range(self.noutputs):
+                for j in range(self.ninputs):
+                    out[i][j] = (polyval(self.num[i][j], x_arr) /
+                                 polyval(self.den[i][j], x_arr))
         return out
 
     def _truncatecoeff(self):
@@ -1048,41 +1055,27 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
         numd, dend, _ = cont2discrete(sys, Twarp, method, alpha)
         return TransferFunction(numd[0, :], dend, Ts)
 
-    def dcgain(self):
+    def dcgain(self, warn_infinite=False):
         """Return the zero-frequency (or DC) gain
 
         For a continous-time transfer function G(s), the DC gain is G(0)
         For a discrete-time transfer function G(z), the DC gain is G(1)
 
+        Parameters
+        ----------
+        warn_infinite : bool, optional
+            By default, don't issue a warning message if the zero-frequency
+            gain is infinite.  Setting `warn_infinite` to generate the warning
+            message.
+
         Returns
         -------
         gain : ndarray
             The zero-frequency gain
-        """
-        if self.isctime():
-            return self._dcgain_cont()
-        else:
-            return self(1)
 
-    def _dcgain_cont(self):
-        """_dcgain_cont() -> DC gain as matrix or scalar
-
-        Special cased evaluation at 0 for continuous-time systems."""
-        gain = np.empty((self.noutputs, self.ninputs), dtype=float)
-        for i in range(self.noutputs):
-            for j in range(self.ninputs):
-                num = self.num[i][j][-1]
-                den = self.den[i][j][-1]
-                if den:
-                    gain[i][j] = num / den
-                else:
-                    if num:
-                        # numerator nonzero: infinite gain
-                        gain[i][j] = np.inf
-                    else:
-                        # numerator is zero too: give up
-                        gain[i][j] = np.nan
-        return np.squeeze(gain)
+        """
+        return self(0, warn_infinite=warn_infinite) if self.isctime() \
+            else self(1, warn_infinite=warn_infinite)
 
     def _isstatic(self):
          """returns True if and only if all of the numerator and denominator