@@ -612,11 +612,14 @@ def dcgain(self):
             at the origin
         """
         try:
-            gain = np.asarray(self.D -
-                              self.C.dot(np.linalg.solve(self.A, self.B)))
+            if self.isctime():
+                gain = np.asarray(self.D -
+                                    self.C.dot(np.linalg.solve(self.A, self.B)))
+            else:
+                gain = self.horner(1)
         except LinAlgError:
-            # zero eigenvalue: singular matrix
-            return np.nan
+            # eigenvalue at DC
+            gain = np.tile(np.nan,(self.outputs,self.inputs))
         return np.squeeze(gain)
 
 
 
@@ -228,7 +228,8 @@ def testArrayAccessSS(self):
 
         assert sys1.dt == sys1_11.dt
 
-    def test_dcgain(self):
+    def test_dcgain_cont(self):
+        """Test DC gain for continuous-time state-space systems"""
         sys = StateSpace(-2.,6.,5.,0)
         np.testing.assert_equal(sys.dcgain(), 15.)
 
@@ -239,6 +240,43 @@ def test_dcgain(self):
         sys3 = StateSpace(0., 1., 1., 0.)
         np.testing.assert_equal(sys3.dcgain(), np.nan)
 
+    def test_dcgain_discr(self):
+        """Test DC gain for discrete-time state-space systems"""
+        # static gain
+        sys = StateSpace([], [], [], 2, True)
+        np.testing.assert_equal(sys.dcgain(), 2)
+
+        # averaging filter
+        sys = StateSpace(0.5, 0.5, 1, 0, True)
+        np.testing.assert_almost_equal(sys.dcgain(), 1)
+
+        # differencer
+        sys = StateSpace(0, 1, -1, 1, True)
+        np.testing.assert_equal(sys.dcgain(), 0)
+
+        # summer
+        sys = StateSpace(1, 1, 1, 0, True)
+        np.testing.assert_equal(sys.dcgain(), np.nan)
+
+    def test_dcgain_integrator(self):
+        """DC gain when eigenvalue at DC returns appropriately sized array of nan"""
+        # the SISO case is also tested in test_dc_gain_{cont,discr}
+        import itertools
+        # iterate over input and output sizes, and continuous (dt=None) and discrete (dt=True) time
+        for inputs,outputs,dt in itertools.product(range(1,6),range(1,6),[None,True]):
+            states = max(inputs,outputs)
+
+            # a matrix that is singular at DC, and has no "useless" states as in _remove_useless_states
+            a = np.triu(np.tile(2,(states,states)))
+            # eigenvalues all +2, except for ...
+            a[0,0] = 0 if dt is None else 1
+            b = np.eye(max(inputs,states))[:states,:inputs]
+            c = np.eye(max(outputs,states))[:outputs,:states]
+            d = np.zeros((outputs,inputs))
+            sys = StateSpace(a,b,c,d,dt)
+            dc = np.squeeze(np.tile(np.nan,(outputs,inputs)))
+            np.testing.assert_array_equal(dc, sys.dcgain())
+
 
     def test_scalarStaticGain(self):
         """Regression: can we create a scalar static gain?"""
 
@@ -524,7 +524,8 @@ def testMatrixMult(self):
         np.testing.assert_array_almost_equal(H.num[1][0], H2.num[0][0])
         np.testing.assert_array_almost_equal(H.den[1][0], H2.den[0][0])
 
-    def test_dcgain(self):
+    def test_dcgain_cont(self):
+        """Test DC gain for continuous-time transfer functions"""
         sys = TransferFunction(6, 3)
         np.testing.assert_equal(sys.dcgain(), 2)
 
@@ -540,6 +541,26 @@ def test_dcgain(self):
         expected = [[5, 7, 11], [2, 2, 2]]
         np.testing.assert_array_equal(sys4.dcgain(), expected)
 
+    def test_dcgain_discr(self):
+        """Test DC gain for discrete-time transfer functions"""
+        # static gain
+        sys = TransferFunction(6, 3, True)
+        np.testing.assert_equal(sys.dcgain(), 2)
+
+        # averaging filter
+        sys = TransferFunction(0.5, [1, -0.5], True)
+        np.testing.assert_almost_equal(sys.dcgain(), 1)
+
+        # differencer
+        sys = TransferFunction(1, [1, -1], True)
+        np.testing.assert_equal(sys.dcgain(), np.inf)
+
+        # summer
+        # causes a RuntimeWarning due to the divide by zero
+        sys = TransferFunction([1,-1], [1], True)
+        np.testing.assert_equal(sys.dcgain(), 0)
+
+
 def suite():
     return unittest.TestLoader().loadTestsFromTestCase(TestXferFcn)
 
 
@@ -945,13 +945,23 @@ def sample(self, Ts, method='zoh', alpha=None):
     def dcgain(self):
         """Return the zero-frequency (or DC) gain
 
-        For a transfer function G(s), the DC gain is G(0)
+        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)
 
         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.outputs, self.inputs), dtype=float)
         for i in range(self.outputs):
             for j in range(self.inputs):