diff --git a/control/dtime.py b/control/dtime.py
index 2ae482811..9b91eabd3 100644
--- a/control/dtime.py
+++ b/control/dtime.py
@@ -55,8 +55,7 @@
 # Sample a continuous time system
 def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None,
         name=None, copy_names=True, **kwargs):
-    """
-    Convert a continuous time system to discrete time by sampling.
+    """Convert a continuous time system to discrete time by sampling.
 
     Parameters
     ----------
@@ -67,9 +66,9 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None,
     method : string
         Method to use for conversion, e.g. 'bilinear', 'zoh' (default)
     alpha : float within [0, 1]
-            The generalized bilinear transformation weighting parameter, which
-            should only be specified with method="gbt", and is ignored
-            otherwise. See :func:`scipy.signal.cont2discrete`.
+        The generalized bilinear transformation weighting parameter, which
+        should only be specified with method="gbt", and is ignored
+        otherwise. See :func:`scipy.signal.cont2discrete`.
     prewarp_frequency : float within [0, infinity)
         The frequency [rad/s] at which to match with the input continuous-
         time system's magnitude and phase (only valid for method='bilinear',
diff --git a/control/freqplot.py b/control/freqplot.py
index 164ec28a2..533515415 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -894,8 +894,10 @@ def gen_zero_centered_series(val_min, val_max, period):
     # list of systems (e.g., "Step response for sys[1], sys[2]").
     #
 
-    # Set the initial title for the data (unique system names)
-    sysnames = list(set([response.sysname for response in data]))
+    # Set the initial title for the data (unique system names, preserving order)
+    seen = set()
+    sysnames = [response.sysname for response in data \
+                if not (response.sysname in seen or seen.add(response.sysname))]
     if title is None:
         if data[0].title is None:
             title = "Bode plot for " + ", ".join(sysnames)
diff --git a/control/tests/discrete_test.py b/control/tests/discrete_test.py
index 96777011e..cccb53708 100644
--- a/control/tests/discrete_test.py
+++ b/control/tests/discrete_test.py
@@ -5,11 +5,12 @@
 
 import numpy as np
 import pytest
+import cmath
 
-from control import (StateSpace, TransferFunction, bode, common_timebase,
-                     feedback, forced_response, impulse_response,
-                     isctime, isdtime, rss, c2d, sample_system, step_response,
-                     timebase)
+import control as ct
+from control import StateSpace, TransferFunction, bode, common_timebase, \
+    feedback, forced_response, impulse_response, isctime, isdtime, rss, \
+    c2d, sample_system, step_response, timebase
 
 
 class TestDiscrete:
@@ -526,3 +527,33 @@ def test_signal_names(self, tsys):
         assert sysd_newnames.find_input('u') is None
         assert sysd_newnames.find_output('y') == 0
         assert sysd_newnames.find_output('x') is None
+
+
+@pytest.mark.parametrize("num, den", [
+    ([1], [1, 1]),
+    ([1, 2], [1, 3]),
+    ([1, 2], [3, 4, 5])
+])
+@pytest.mark.parametrize("dt", [True, 0.1, 2])
+@pytest.mark.parametrize("method", ['zoh', 'bilinear', 'matched'])
+def test_c2d_matched(num, den, dt, method):
+    sys_ct = ct.tf(num, den)
+    sys_dt = ct.sample_system(sys_ct, dt, method=method)
+    assert sys_dt.dt == dt                      # make sure sampling time is OK
+    assert cmath.isclose(sys_ct(0), sys_dt(1))  # check zero frequency gain
+    assert cmath.isclose(
+        sys_ct.dcgain(), sys_dt.dcgain())       # another way to check
+
+    if method in ['zoh', 'matched']:
+        # Make sure that poles were properly matched
+        zpoles = sys_dt.poles()
+        for cpole in sys_ct.poles():
+            zpole = zpoles[(np.abs(zpoles - cmath.exp(cpole * dt))).argmin()]
+            assert cmath.isclose(cmath.exp(cpole * dt), zpole)
+
+    if method in ['matched']:
+        # Make sure that zeros were properly matched
+        zzeros = sys_dt.zeros()
+        for czero in sys_ct.zeros():
+            zzero = zzeros[(np.abs(zzeros - cmath.exp(czero * dt))).argmin()]
+            assert cmath.isclose(cmath.exp(czero * dt), zzero)
diff --git a/control/xferfcn.py b/control/xferfcn.py
index b5334b7b8..206f0bc92 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1181,7 +1181,9 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None,
         if not self.issiso():
             raise ControlMIMONotImplemented("Not implemented for MIMO systems")
         if method == "matched":
-            return _c2d_matched(self, Ts)
+            if prewarp_frequency is not None:
+                warn('prewarp_frequency ignored: incompatible conversion')
+            return _c2d_matched(self, Ts, name=name, **kwargs)
         sys = (self.num[0][0], self.den[0][0])
         if prewarp_frequency is not None:
             if method in ('bilinear', 'tustin') or \
@@ -1284,9 +1286,12 @@ def _isstatic(self):
 
 
 # c2d function contributed by Benjamin White, Oct 2012
-def _c2d_matched(sysC, Ts):
+def _c2d_matched(sysC, Ts, **kwargs):
+    if not sysC.issiso():
+        raise ControlMIMONotImplemented("Not implemented for MIMO systems")
+
     # Pole-zero match method of continuous to discrete time conversion
-    szeros, spoles, sgain = tf2zpk(sysC.num[0][0], sysC.den[0][0])
+    szeros, spoles, _ = tf2zpk(sysC.num[0][0], sysC.den[0][0])
     zzeros = [0] * len(szeros)
     zpoles = [0] * len(spoles)
     pregainnum = [0] * len(szeros)
@@ -1302,9 +1307,9 @@ def _c2d_matched(sysC, Ts):
         zpoles[idx] = z
         pregainden[idx] = 1 - z
     zgain = np.multiply.reduce(pregainnum) / np.multiply.reduce(pregainden)
-    gain = sgain / zgain
+    gain = sysC.dcgain() / zgain
     sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain)
-    return TransferFunction(sysDnum, sysDden, Ts)
+    return TransferFunction(sysDnum, sysDden, Ts, **kwargs)
 
 
 # Utility function to convert a transfer function polynomial to a string