diff --git a/control/bench/time_freqresp.py b/control/bench/time_freqresp.py
index 1945cbc24..3ae837082 100644
--- a/control/bench/time_freqresp.py
+++ b/control/bench/time_freqresp.py
@@ -8,7 +8,7 @@
 sys_tf = tf(sys)
 w = logspace(-1,1,50)
 ntimes = 1000
-time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes)
-time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes)
+time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes)
+time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes)
 print("State-space model on %d states: %f" % (nstates, time_ss))
 print("Transfer-function model on %d states: %f" % (nstates, time_tf))
diff --git a/control/frdata.py b/control/frdata.py
index 8ca9dfd9d..95c84e0ba 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -48,7 +48,7 @@
 from warnings import warn
 import numpy as np
 from numpy import angle, array, empty, ones, \
-    real, imag, absolute, eye, linalg, where, dot
+    real, imag, absolute, eye, linalg, where, dot, sort
 from scipy.interpolate import splprep, splev
 from .lti import LTI
 
@@ -100,6 +100,7 @@ def __init__(self, *args, **kwargs):
         object, other than an FRD, call FRD(sys, omega)
 
         """
+        # TODO: discrete-time FRD systems?
         smooth = kwargs.get('smooth', False)
 
         if len(args) == 2:
@@ -107,17 +108,10 @@ def __init__(self, *args, **kwargs):
                 # not an FRD, but still a system, second argument should be
                 # the frequency range
                 otherlti = args[0]
-                self.omega = array(args[1], dtype=float)
-                self.omega.sort()
+                self.omega = sort(np.asarray(args[1], dtype=float))
                 numfreq = len(self.omega)
-
                 # calculate frequency response at my points
-                self.fresp = empty(
-                    (otherlti.outputs, otherlti.inputs, numfreq),
-                    dtype=complex)
-                for k, w in enumerate(self.omega):
-                    self.fresp[:, :, k] = otherlti._evalfr(w)
-
+                self.fresp = otherlti(1j * self.omega, squeeze=False)
             else:
                 # The user provided a response and a freq vector
                 self.fresp = array(args[0], dtype=complex)
@@ -141,7 +135,7 @@ def __init__(self, *args, **kwargs):
             self.omega = args[0].omega
             self.fresp = args[0].fresp
         else:
-            raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
+            raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
         # create interpolation functions
         if smooth:
@@ -163,7 +157,6 @@ def __str__(self):
         mimo = self.inputs > 1 or self.outputs > 1
         outstr = ['Frequency response data']
 
-        mt, pt, wt = self.freqresp(self.omega)
         for i in range(self.inputs):
             for j in range(self.outputs):
                 if mimo:
@@ -171,9 +164,10 @@ def __str__(self):
                 outstr.append('Freq [rad/s]  Response')
                 outstr.append('------------  ---------------------')
                 outstr.extend(
-                    ['%12.3f  %10.4g%+10.4gj' % (w, m, p)
-                     for m, p, w in zip(real(self.fresp[j, i, :]),
-                                        imag(self.fresp[j, i, :]), wt)])
+                    ['%12.3f  %10.4g%+10.4gj' % (w, re, im)
+                     for w, re, im in zip(self.omega,
+                                          real(self.fresp[j, i, :]),
+                                          imag(self.fresp[j, i, :]))])
 
         return '\n'.join(outstr)
 
@@ -342,110 +336,115 @@ def __pow__(self, other):
             return (FRD(ones(self.fresp.shape), self.omega) / self) * \
                 (self**(other+1))
 
-    def evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self._evalfr(omega) returns the value of the frequency response
-        at frequency omega.
-
-        Note that a "normal" FRD only returns values for which there is an
-        entry in the omega vector. An interpolating FRD can return
-        intermediate values.
-
-        """
-        warn("FRD.evalfr(omega) will be deprecated in a future release "
-             "of python-control; use sys.eval(omega) instead",
-             PendingDeprecationWarning)         # pragma: no coverage
-        return self._evalfr(omega)
-
     # Define the `eval` function to evaluate an FRD at a given (real)
     # frequency.  Note that we choose to use `eval` instead of `evalfr` to
     # avoid confusion with :func:`evalfr`, which takes a complex number as its
     # argument.  Similarly, we don't use `__call__` to avoid confusion between
     # G(s) for a transfer function and G(omega) for an FRD object.
-    def eval(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self.evalfr(omega) returns the value of the frequency response
-        at frequency omega.
+    # update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform
+    # interface to systems in general and the lti.frequency_response method
+    def eval(self, omega, squeeze=True):
+        """Evaluate a transfer function at angular frequency omega.
 
         Note that a "normal" FRD only returns values for which there is an
         entry in the omega vector. An interpolating FRD can return
         intermediate values.
 
+        Parameters
+        ----------
+        omega: float or array_like
+            Frequencies in radians per second
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of `omega`
+
+        Returns
+        -------
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is ``len(x)`` if and
+            only if system is SISO and ``squeeze=True``.
+
         """
-        return self._evalfr(omega)
-
-    # Internal function to evaluate the frequency responses
-    def _evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency."""
-        # Preallocate the output.
-        if getattr(omega, '__iter__', False):
-            out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
-        else:
-            out = empty((self.outputs, self.inputs), dtype=complex)
+        omega_array = np.array(omega, ndmin=1) # array-like version of omega
+        if any(omega_array.imag > 0):
+            raise ValueError("FRD.eval can only accept real-valued omega")
 
         if self.ifunc is None:
-            try:
-                out = self.fresp[:, :, where(self.omega == omega)[0][0]]
-            except Exception:
+            elements = np.isin(self.omega, omega) # binary array
+            if sum(elements) < len(omega_array):
                 raise ValueError(
-                    "Frequency %f not in frequency list, try an interpolating"
-                    " FRD if you want additional points" % omega)
-        else:
-            if getattr(omega, '__iter__', False):
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        for k, w in enumerate(omega):
-                            frraw = splev(w, self.ifunc[i, j], der=0)
-                            out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+                    "not all frequencies omega are in frequency list of FRD "
+                    "system. Try an interpolating FRD for additional points.")
             else:
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        frraw = splev(omega, self.ifunc[i, j], der=0)
-                        out[i, j] = frraw[0] + 1.0j * frraw[1]
-
+                out = self.fresp[:, :, elements]
+        else:
+            out = empty((self.outputs, self.inputs, len(omega_array)),
+                         dtype=complex)
+            for i in range(self.outputs):
+                for j in range(self.inputs):
+                    for k, w in enumerate(omega_array):
+                        frraw = splev(w, self.ifunc[i, j], der=0)
+                        out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+        if not hasattr(omega, '__len__'):
+            # omega is a scalar, squeeze down array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            out = out[0][0]
         return out
 
-    # Method for generating the frequency response of the system
-    def freqresp(self, omega):
-        """Evaluate the frequency response at a list of angular frequencies.
+    def __call__(self, s, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies.
 
-        Reports the value of the magnitude, phase, and angular frequency of
-        the requency response evaluated at omega, where omega is a list of
-        angular frequencies, and is a sorted version of the input omega.
+        Returns the complex frequency response `sys(s)`.
+
+        To evaluate at a frequency omega in radians per second, enter
+        ``x = omega * 1j`` or use ``sys.eval(omega)``
 
         Parameters
         ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
+        s: complex scalar or array_like
+            Complex frequencies
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of x`.
 
         Returns
         -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray or list or tuple
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # Preallocate outputs.
-        numfreq = len(omega)
-        mag = empty((self.outputs, self.inputs, numfreq))
-        phase = empty((self.outputs, self.inputs, numfreq))
+        fresp : (self.outputs, self.inputs, len(s)) or len(s) complex ndarray
+            The frequency response of the system. Array is ``len(s)`` if and
+            only if system is SISO and ``squeeze=True``.
 
-        omega.sort()
+        Raises
+        ------
+        ValueError
+            If `s` is not purely imaginary, because
+            :class:`FrequencyDomainData` systems are only defined at imaginary
+            frequency values.
+
+        """
+        if any(abs(np.array(s, ndmin=1).real) > 0):
+            raise ValueError("__call__: FRD systems can only accept"
+                            "purely imaginary frequencies")
+        # need to preserve array or scalar status
+        if hasattr(s, '__len__'):
+            return self.eval(np.asarray(s).imag, squeeze=squeeze)
+        else:
+            return self.eval(complex(s).imag, squeeze=squeeze)
 
-        for k, w in enumerate(omega):
-            fresp = self._evalfr(w)
-            mag[:, :, k] = abs(fresp)
-            phase[:, :, k] = angle(fresp)
+    def freqresp(self, omega):
+        """(deprecated) Evaluate transfer function at complex frequencies.
 
-        return mag, phase, omega
+        .. deprecated::0.9.0
+            Method has been given the more pythonic name
+            :meth:`FrequencyResponseData.frequency_response`. Or use
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("FrequencyResponseData.freqresp(omega) will be removed in a "
+             "future release of python-control; use "
+             "FrequencyResponseData.frequency_response(omega), or "
+             "freqresp(sys, omega) in the MATLAB compatibility module "
+             "instead", DeprecationWarning)
+        return self.frequency_response(omega)
 
     def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two FRD objects."""
@@ -515,11 +514,10 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
             "Frequency ranges of FRD do not match, conversion not implemented")
 
     elif isinstance(sys, LTI):
-        omega.sort()
-        fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
-        for k, w in enumerate(omega):
-            fresp[:, :, k] = sys._evalfr(w)
-
+        omega = np.sort(omega)
+        fresp = sys(1j * omega)
+        if len(fresp.shape) == 1:
+            fresp = fresp[np.newaxis, np.newaxis, :]
         return FRD(fresp, omega, smooth=True)
 
     elif isinstance(sys, (int, float, complex, np.number)):
diff --git a/control/freqplot.py b/control/freqplot.py
index 5a03694db..38b525aec 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -151,14 +151,14 @@ def bode_plot(syslist, omega=None,
 
     Notes
     -----
-    1. Alternatively, you may use the lower-level method
-       ``(mag, phase, freq) = sys.freqresp(freq)`` to generate the frequency
-       response for a system, but it returns a MIMO response.
+    1. Alternatively, you may use the lower-level methods
+       :meth:`LTI.frequency_response` or ``sys(s)`` or ``sys(z)`` or to
+       generate the frequency response for a single system.
 
     2. If a discrete time model is given, the frequency response is plotted
-       along the upper branch of the unit circle, using the mapping z = exp(j
-       \\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
-       timebase.  If not timebase is specified (dt = True), dt is set to 1.
+       along the upper branch of the unit circle, using the mapping ``z = exp(1j
+       * omega * dt)`` where `omega` ranges from 0 to `pi/dt` and `dt` is the discrete
+       timebase.  If timebase not specified (``dt=True``), `dt` is set to 1.
 
     Examples
     --------
@@ -228,7 +228,7 @@ def bode_plot(syslist, omega=None,
                 nyquistfrq = None
 
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
+            mag_tmp, phase_tmp, omega_sys = sys.frequency_response(omega_sys)
             mag = np.atleast_1d(np.squeeze(mag_tmp))
             phase = np.atleast_1d(np.squeeze(phase_tmp))
 
@@ -588,7 +588,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
                 "Nyquist is currently only implemented for SISO systems.")
         else:
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+            mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
             mag = np.squeeze(mag_tmp)
             phase = np.squeeze(phase_tmp)
 
@@ -718,7 +718,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     omega_plot = omega / (2. * math.pi) if Hz else omega
 
     # TODO: Need to add in the mag = 1 lines
-    mag_tmp, phase_tmp, omega = S.freqresp(omega)
+    mag_tmp, phase_tmp, omega = S.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['s'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -728,7 +728,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['s'].tick_params(labelbottom=False)
     plot_axes['s'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (P * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (P * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['ps'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -738,7 +738,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['ps'].set_ylabel("$|PS|$" + " (dB)" if dB else "")
     plot_axes['ps'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (C * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (C * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['cs'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -749,7 +749,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['cs'].set_ylabel("$|CS|$" + " (dB)" if dB else "")
     plot_axes['cs'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = T.freqresp(omega)
+    mag_tmp, phase_tmp, omega = T.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['t'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
diff --git a/control/lti.py b/control/lti.py
index e41fe416b..eff14be2a 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -13,11 +13,11 @@
 """
 
 import numpy as np
-from numpy import absolute, real
+from numpy import absolute, real, angle, abs
 from warnings import warn
 
-__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual', 
-           'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr', 
+__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual',
+           'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr',
            'freqresp', 'dcgain']
 
 class LTI:
@@ -111,11 +111,56 @@ def damp(self):
         Z = -real(splane_poles)/wn
         return wn, Z, poles
 
+    def frequency_response(self, omega, squeeze=True):
+        """Evaluate the linear time-invariant system at an array of angular
+        frequencies.
+
+        Reports the frequency response of the system,
+
+             G(j*omega) = mag*exp(j*phase)
+
+        for continuous time systems. For discrete time systems, the response is
+        evaluated around the unit circle such that
+
+             G(exp(j*omega*dt)) = mag*exp(j*phase).
+
+        Parameters
+        ----------
+        omega : array_like or float
+            A list, tuple, array, or scalar value of frequencies in
+            radians/sec at which the system will be evaluated.
+        squeeze: bool, optional (default=True)
+            If True and sys is single input, single output (SISO), return a
+            1D array or scalar depending on omega's length.
+
+        Returns
+        -------
+        mag : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
+            The magnitude (absolute value, not dB or log10) of the system
+            frequency response.
+        phase : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
+            The wrapped phase in radians of the system frequency response.
+        omega : ndarray
+            The (sorted) frequencies at which the response was evaluated.
+
+        """
+        omega = np.sort(np.array(omega, ndmin=1))
+        if isdtime(self, strict=True):
+            # Convert the frequency to discrete time
+            if np.any(omega * self.dt > np.pi):
+                warn("__call__: evaluation above Nyquist frequency")
+            s = np.exp(1j * omega * self.dt)
+        else:
+            s = 1j * omega
+        response = self.__call__(s, squeeze=squeeze)
+        return abs(response), angle(response), omega
+
     def dcgain(self):
         """Return the zero-frequency gain"""
         raise NotImplementedError("dcgain not implemented for %s objects" %
                                   str(self.__class__))
 
+
 # Test to see if a system is SISO
 def issiso(sys, strict=False):
     """
@@ -162,50 +207,50 @@ def timebase(sys, strict=True):
 def common_timebase(dt1, dt2):
     """
     Find the common timebase when interconnecting systems
-    
+
     Parameters
     ----------
-    dt1, dt2: number or system with a 'dt' attribute (e.g. TransferFunction 
+    dt1, dt2: number or system with a 'dt' attribute (e.g. TransferFunction
         or StateSpace system)
 
     Returns
     -------
     dt: number
-        The common timebase of dt1 and dt2, as specified in 
-        :ref:`conventions-ref`. 
-        
+        The common timebase of dt1 and dt2, as specified in
+        :ref:`conventions-ref`.
+
     Raises
     ------
     ValueError
         when no compatible time base can be found
     """
-    # explanation: 
+    # explanation:
     # if either dt is None, they are compatible with anything
-    # if either dt is True (discrete with unspecified time base), 
+    # if either dt is True (discrete with unspecified time base),
     #   use the timebase of the other, if it is also discrete
-    # otherwise both dts must be equal 
+    # otherwise both dts must be equal
     if hasattr(dt1, 'dt'):
         dt1 = dt1.dt
     if hasattr(dt2, 'dt'):
         dt2 = dt2.dt
 
-    if dt1 is None: 
+    if dt1 is None:
         return dt2
-    elif dt2 is None: 
+    elif dt2 is None:
         return dt1
-    elif dt1 is True: 
+    elif dt1 is True:
         if dt2 > 0:
             return dt2
-        else: 
+        else:
             raise ValueError("Systems have incompatible timebases")
-    elif dt2 is True: 
-        if dt1 > 0: 
+    elif dt2 is True:
+        if dt1 > 0:
             return dt1
-        else: 
+        else:
             raise ValueError("Systems have incompatible timebases")
     elif np.isclose(dt1, dt2):
         return dt1
-    else: 
+    else:
         raise ValueError("Systems have incompatible timebases")
 
 # Check to see if two timebases are equal
@@ -221,9 +266,9 @@ def timebaseEqual(sys1, sys2):
     timebase (dt > 0) then their timebases must be equal.
     """
     warn("timebaseEqual will be deprecated in a future release of "
-         "python-control; use :func:`common_timebase` instead", 
+         "python-control; use :func:`common_timebase` instead",
          PendingDeprecationWarning)
-        
+
     if (type(sys1.dt) == bool or type(sys2.dt) == bool):
         # Make sure both are unspecified discrete timebases
         return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt
@@ -413,24 +458,33 @@ def damp(sys, doprint=True):
                       (p.real, p.imag, d, w))
     return wn, damping, poles
 
-def evalfr(sys, x):
+def evalfr(sys, x, squeeze=True):
     """
-    Evaluate the transfer function of an LTI system for a single complex
-    number x.
+    Evaluate the transfer function of an LTI system for complex frequency x.
+
+    Returns the complex frequency response `sys(x)` where `x` is `s` for
+    continuous-time systems and `z` for discrete-time systems.
 
-    To evaluate at a frequency, enter x = omega*j, where omega is the
-    frequency in radians
+    To evaluate at a frequency omega in radians per second, enter
+    ``x = omega * 1j`` for continuous-time systems, or
+    ``x = exp(1j * omega * dt)`` for discrete-time systems, or use
+    ``freqresp(sys, omega)``.
 
     Parameters
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    x: scalar
-        Complex number
+    x: complex scalar or array_like
+        Complex frequency(s)
+    squeeze: bool, optional (default=True)
+        If True and sys is single input single output (SISO), returns a
+        1D array or scalar depending on the length of x.
 
     Returns
     -------
-    fresp: ndarray
+    fresp : (sys.outputs, sys.inputs, len(x)) or len(x) complex ndarray
+        The frequency response of the system. Array is len(x) if and only if
+        system is SISO and squeeze=True.
 
     See Also
     --------
@@ -439,8 +493,8 @@ def evalfr(sys, x):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.evalfr and
-    TransferFunction.evalfr.
+    This function is a wrapper for StateSpace.__call__ and
+    TransferFunction.__call__.
 
     Examples
     --------
@@ -451,12 +505,9 @@ def evalfr(sys, x):
 
     .. todo:: Add example with MIMO system
     """
-    if issiso(sys):
-        return sys.horner(x)[0][0]
-    return sys.horner(x)
-
+    return sys.__call__(x, squeeze=squeeze)
 
-def freqresp(sys, omega):
+def freqresp(sys, omega, squeeze=True):
     """
     Frequency response of an LTI system at multiple angular frequencies.
 
@@ -464,19 +515,22 @@ def freqresp(sys, omega):
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    omega: array_like
+    omega: float or array_like
         A list of frequencies in radians/sec at which the system should be
         evaluated. The list can be either a python list or a numpy array
         and will be sorted before evaluation.
+    squeeze: bool, optional (default=True)
+        If True and sys is single input, single output (SISO), returns
+        1D array or scalar depending on omega's length.
 
     Returns
     -------
-    mag : (self.outputs, self.inputs, len(omega)) ndarray
+    mag : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray
         The magnitude (absolute value, not dB or log10) of the system
         frequency response.
-    phase : (self.outputs, self.inputs, len(omega)) ndarray
+    phase : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray
         The wrapped phase in radians of the system frequency response.
-    omega : ndarray or list or tuple
+    omega : ndarray
         The list of sorted frequencies at which the response was
         evaluated.
 
@@ -487,9 +541,8 @@ def freqresp(sys, omega):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.freqresp and
-    TransferFunction.freqresp.  The output omega is a sorted version of the
-    input omega.
+    This function is a wrapper for StateSpace.frequency_response and
+    TransferFunction.frequency_response.
 
     Examples
     --------
@@ -514,7 +567,7 @@ def freqresp(sys, omega):
         #>>> # frequency response from the 1st input to the 2nd output, for
         #>>> # s = 0.1i, i, 10i.
     """
-    return sys.freqresp(omega)
+    return sys.frequency_response(omega, squeeze=squeeze)
 
 
 def dcgain(sys):
diff --git a/control/margins.py b/control/margins.py
index 03e78352f..20da2a879 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -339,24 +339,24 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
         # a bit coarse, have the interpolated frd evaluated again
         def _mod(w):
             """Calculate |G(jw)| - 1"""
-            return np.abs(sys._evalfr(w)[0][0]) - 1
+            return np.abs(sys(1j * w)) - 1
 
         def _arg(w):
             """Calculate the phase angle at -180 deg"""
-            return np.angle(-sys._evalfr(w)[0][0])
+            return np.angle(-sys(1j * w))
 
         def _dstab(w):
             """Calculate the distance from -1 point"""
-            return np.abs(sys._evalfr(w)[0][0] + 1.)
+            return np.abs(sys(1j * w) + 1.)
 
         # find the phase crossings ang(H(jw) == -180
         widx = np.where(np.diff(np.sign(_arg(sys.omega))))[0]
-        widx = widx[np.real(sys._evalfr(sys.omega[widx])[0][0]) <= 0]
+        widx = widx[np.real(sys(1j * sys.omega[widx])) <= 0]
         w_180 = np.array(
             [sp.optimize.brentq(_arg, sys.omega[i], sys.omega[i+1])
              for i in widx])
         # TODO: replace by evalfr(sys, 1J*w) or sys(1J*w), (needs gh-449)
-        w180_resp = sys._evalfr(w_180)[0][0]
+        w180_resp = sys(1j * w_180)
 
         # Find all crossings, note that this depends on omega having
         # a correct range
@@ -364,7 +364,7 @@ def _dstab(w):
         wc = np.array(
             [sp.optimize.brentq(_mod, sys.omega[i], sys.omega[i+1])
              for i in widx])
-        wc_resp = sys._evalfr(wc)[0][0]
+        wc_resp = sys(1j * wc)
 
         # find all stab margins?
         widx, = np.where(np.diff(np.sign(np.diff(_dstab(sys.omega)))) > 0)
@@ -374,7 +374,7 @@ def _dstab(w):
                                          ).x
              for i in widx])
         wstab = wstab[(wstab >= sys.omega[0]) * (wstab <= sys.omega[-1])]
-        ws_resp = sys._evalfr(wstab)[0][0]
+        ws_resp = sys(1j * wstab)
 
     with np.errstate(all='ignore'):  # |G|=0 is okay and yields inf
         GM = 1. / np.abs(w180_resp)
diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py
index 413dc6d86..6b88214c6 100644
--- a/control/matlab/__init__.py
+++ b/control/matlab/__init__.py
@@ -224,8 +224,8 @@
 \*  :func:`~control.nichols`    Nichols plot
 \*  :func:`margin`              gain and phase margins
 \   lti/allmargin               all crossover frequencies and margins
-\*  :func:`freqresp`            frequency response over a frequency grid
-\*  :func:`evalfr`              frequency response at single frequency
+\*  :func:`freqresp`            frequency response
+\*  :func:`evalfr`              frequency response at complex frequency s
 ==  ==========================  ============================================
 
 
diff --git a/control/nichols.py b/control/nichols.py
index ca0505957..c1d8ff9b6 100644
--- a/control/nichols.py
+++ b/control/nichols.py
@@ -95,7 +95,7 @@ def nichols_plot(sys_list, omega=None, grid=None):
 
     for sys in sys_list:
         # Get the magnitude and phase of the system
-        mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+        mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
         mag = np.squeeze(mag_tmp)
         phase = np.squeeze(phase_tmp)
 
diff --git a/control/rlocus.py b/control/rlocus.py
index 479a833ab..3a3c1c2b6 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -476,7 +476,7 @@ def _systopoly1d(sys):
         sys = _convert_to_transfer_function(sys)
 
         # Make sure we have a SISO system
-        if (sys.inputs > 1 or sys.outputs > 1):
+        if not sys.issiso():
             raise ControlMIMONotImplemented()
 
         # Start by extracting the numerator and denominator from system object
@@ -497,7 +497,7 @@ def _RLFindRoots(nump, denp, kvect):
     """Find the roots for the root locus."""
     # Convert numerator and denominator to polynomials if they aren't
     roots = []
-    for k in kvect:
+    for k in np.array(kvect, ndmin=1):
         curpoly = denp + k * nump
         curroots = curpoly.r
         if len(curroots) < denp.order:
@@ -581,10 +581,10 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
     # Catch type error when event click is in the figure but not in an axis
     try:
         s = complex(event.xdata, event.ydata)
-        K = -1. / sys.horner(s)
-        K_xlim = -1. / sys.horner(
+        K = -1. / sys(s)
+        K_xlim = -1. / sys(
             complex(event.xdata + 0.05 * abs(xlim[1] - xlim[0]), event.ydata))
-        K_ylim = -1. / sys.horner(
+        K_ylim = -1. / sys(
             complex(event.xdata, event.ydata + 0.05 * abs(ylim[1] - ylim[0])))
 
     except TypeError:
@@ -625,7 +625,7 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
             ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8,
                            zorder=20, label='gain_point')
 
-        return K.real[0][0]
+        return K.real
 
 
 def _removeLine(label, ax):
diff --git a/control/statesp.py b/control/statesp.py
index b41f18c7d..8bf4c0d99 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -54,7 +54,7 @@
 import math
 import numpy as np
 from numpy import any, array, asarray, concatenate, cos, delete, \
-    dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze
+    dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze, pi
 from numpy.random import rand, randn
 from numpy.linalg import solve, eigvals, matrix_rank
 from numpy.linalg.linalg import LinAlgError
@@ -640,125 +640,146 @@ def __rdiv__(self, other):
         raise NotImplementedError(
             "StateSpace.__rdiv__ is not implemented yet.")
 
-    def evalfr(self, omega):
-        """Evaluate a SS system's transfer function at a single frequency.
+    def __call__(self, x, squeeze=True):
+        """Evaluate system's transfer function at complex frequency.
 
-        self._evalfr(omega) returns the value of the transfer function matrix
-        with input value s = i * omega.
+        Returns the complex frequency response `sys(x)` where `x` is `s` for
+        continuous-time systems and `z` for discrete-time systems.
 
-        """
-        warn("StateSpace.evalfr(omega) will be deprecated in a future "
-             "release of python-control; use evalfr(sys, omega*1j) instead",
-             PendingDeprecationWarning)
-        return self._evalfr(omega)
-
-    def _evalfr(self, omega):
-        """Evaluate a SS system's transfer function at a single frequency"""
-        # Figure out the point to evaluate the transfer function
-        if isdtime(self, strict=True):
-            s = exp(1.j * omega * self.dt)
-            if omega * self.dt > math.pi:
-                warn("_evalfr: frequency evaluation above Nyquist frequency")
-        else:
-            s = omega * 1.j
+        To evaluate at a frequency omega in radians per second, enter
+        ``x = omega * 1j``, for continuous-time systems, or
+        ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use
+        :meth:`StateSpace.frequency_response`.
 
-        return self.horner(s)
+        Parameters
+        ----------
+        x: complex or complex array_like
+            Complex frequencies
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of `x`.
 
-    def horner(self, s):
-        """Evaluate the systems's transfer function for a complex variable
+        Returns
+        -------
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is ``len(x)`` if and
+            only if system is SISO and ``squeeze=True``.
 
-        Returns a matrix of values evaluated at complex variable s.
         """
-        resp = np.dot(self.C, solve(s * eye(self.states) - self.A,
-                                    self.B)) + self.D
-        return array(resp)
+        # Use Slycot if available
+        out = self.horner(x)
+        if not hasattr(x, '__len__'):
+            # received a scalar x, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            return out[0][0]
+        else:
+            return out
 
-    def freqresp(self, omega):
-        """Evaluate the system's transfer function at a list of frequencies
+    def slycot_laub(self, x):
+        """Evaluate system's transfer function at complex frequency
+        using Laub's method from Slycot.
+
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
+        for a more user-friendly interface.
+
+        Parameters
+        ----------
+        x : complex array_like or complex
+            Complex frequency
 
-        Reports the frequency response of the system,
+        Returns
+        -------
+        output : (number_outputs, number_inputs, len(x)) complex ndarray
+            Frequency response
+        """
+        from slycot import tb05ad
+
+        # preallocate
+        x_arr = np.atleast_1d(x) # array-like version of x
+        n = self.states
+        m = self.inputs
+        p = self.outputs
+        out = np.empty((p, m, len(x_arr)), dtype=complex)
+        # The first call both evaluates C(sI-A)^-1 B and also returns
+        # Hessenberg transformed matrices at, bt, ct.
+        result = tb05ad(n, m, p, x_arr[0], self.A, self.B, self.C, job='NG')
+        # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
+        at = result[0]
+        bt = result[1]
+        ct = result[2]
+
+        # TB05AD frequency evaluation does not include direct feedthrough.
+        out[:, :, 0] = result[3] + self.D
+
+        # Now, iterate through the remaining frequencies using the
+        # transformed state matrices, at, bt, ct.
+
+        # Start at the second frequency, already have the first.
+        for kk, x_kk in enumerate(x_arr[1:len(x_arr)]):
+            result = tb05ad(n, m, p, x_kk, at, bt, ct, job='NH')
+            # When job='NH', result = (g_i, hinvb, info)
+
+            # kk+1 because enumerate starts at kk = 0.
+            # but zero-th spot is already filled.
+            out[:, :, kk+1] = result[0] + self.D
+        return out
 
-             G(j*omega) = mag*exp(j*phase)
+    def horner(self, x):
+        """Evaluate system's transfer function at complex frequency
+        using Laub's or Horner's method.
 
-        for continuous time. For discrete time systems, the response is
-        evaluated around the unit circle such that
+        Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z`
+        for discrete-time systems.
 
-             G(exp(j*omega*dt)) = mag*exp(j*phase).
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
+        for a more user-friendly interface.
 
         Parameters
         ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
+        x : complex array_like or complex
+            Complex frequencies
 
         Returns
         -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # In case omega is passed in as a list, rather than a proper array.
-        omega = np.asarray(omega)
-
-        numFreqs = len(omega)
-        Gfrf = np.empty((self.outputs, self.inputs, numFreqs),
-                        dtype=np.complex128)
-
-        # Sort frequency and calculate complex frequencies on either imaginary
-        # axis (continuous time) or unit circle (discrete time).
-        omega.sort()
-        if isdtime(self, strict=True):
-            cmplx_freqs = exp(1.j * omega * self.dt)
-            if max(np.abs(omega)) * self.dt > math.pi:
-                warn("freqresp: frequency evaluation above Nyquist frequency")
-        else:
-            cmplx_freqs = omega * 1.j
+        output : (self.outputs, self.inputs, len(x)) complex ndarray
+            Frequency response
 
-        # Do the frequency response evaluation. Use TB05AD from Slycot
-        # if it's available, otherwise use the built-in horners function.
+        Notes
+        -----
+        Attempts to use Laub's method from Slycot library, with a
+        fall-back to python code.
+        """
         try:
-            from slycot import tb05ad
-
-            n = np.shape(self.A)[0]
-            m = self.inputs
-            p = self.outputs
-            # The first call both evaluates C(sI-A)^-1 B and also returns
-            # Hessenberg transformed matrices at, bt, ct.
-            result = tb05ad(n, m, p, cmplx_freqs[0], self.A,
-                            self.B, self.C, job='NG')
-            # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
-            at = result[0]
-            bt = result[1]
-            ct = result[2]
-
-            # TB05AD frequency evaluation does not include direct feedthrough.
-            Gfrf[:, :, 0] = result[3] + self.D
-
-            # Now, iterate through the remaining frequencies using the
-            # transformed state matrices, at, bt, ct.
-
-            # Start at the second frequency, already have the first.
-            for kk, cmplx_freqs_kk in enumerate(cmplx_freqs[1:numFreqs]):
-                result = tb05ad(n, m, p, cmplx_freqs_kk, at,
-                                bt, ct, job='NH')
-                # When job='NH', result = (g_i, hinvb, info)
-
-                # kk+1 because enumerate starts at kk = 0.
-                # but zero-th spot is already filled.
-                Gfrf[:, :, kk+1] = result[0] + self.D
-
-        except ImportError:  # Slycot unavailable. Fall back to horner.
-            for kk, cmplx_freqs_kk in enumerate(cmplx_freqs):
-                Gfrf[:, :, kk] = self.horner(cmplx_freqs_kk)
-
-        #      mag           phase           omega
-        return np.abs(Gfrf), np.angle(Gfrf), omega
+            out = self.slycot_laub(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
+            # Preallocate
+            out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
+
+            #TODO: can this be vectorized?
+            for idx, x_idx in enumerate(x_arr):
+                out[:,:,idx] = \
+                    np.dot(self.C,
+                        solve(x_idx * eye(self.states) - self.A, self.B)) \
+                    + self.D
+        return out
+
+    def freqresp(self, omega):
+        """(deprecated) Evaluate transfer function at complex frequencies.
+
+        .. deprecated::0.9.0
+            Method has been given the more pythonic name
+            :meth:`StateSpace.frequency_response`. Or use
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("StateSpace.freqresp(omega) will be removed in a "
+             "future release of python-control; use "
+             "sys.frequency_response(omega), or freqresp(sys, omega) in the "
+             "MATLAB compatibility module instead", DeprecationWarning)
+        return self.frequency_response(omega)
 
     # Compute poles and zeros
     def pole(self):
@@ -1148,7 +1169,7 @@ def dcgain(self):
                 gain = np.asarray(self.D -
                                   self.C.dot(np.linalg.solve(self.A, self.B)))
             else:
-                gain = self.horner(1)
+                gain = np.squeeze(self.horner(1))
         except LinAlgError:
             # eigenvalue at DC
             gain = np.tile(np.nan, (self.outputs, self.inputs))
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index 18f2f17b1..5343112fe 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -39,8 +39,8 @@ def testSISOtf(self):
         frd = FRD(h, omega)
         assert isinstance(frd, FRD)
 
-        mag1, phase1, omega1 = frd.freqresp([1.0])
-        mag2, phase2, omega2 = h.freqresp([1.0])
+        mag1, phase1, omega1 = frd.frequency_response([1.0])
+        mag2, phase2, omega2 = h.frequency_response([1.0])
         np.testing.assert_array_almost_equal(mag1, mag2)
         np.testing.assert_array_almost_equal(phase1, phase2)
         np.testing.assert_array_almost_equal(omega1, omega2)
@@ -54,39 +54,39 @@ def testOperators(self):
         f2 = FRD(h2, omega)
 
         np.testing.assert_array_almost_equal(
-            (f1 + f2).freqresp([0.1, 1.0, 10])[0],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 + f2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 + f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 + f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 - f2).freqresp([0.1, 1.0, 10])[0],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 - f2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 - f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 - f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[1])
 
         # multiplication and division
         np.testing.assert_array_almost_equal(
-            (f1 * f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 * h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 * f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 / h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 / f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / h2).frequency_response([0.1, 1.0, 10])[1])
 
         # with default conversion from scalar
         np.testing.assert_array_almost_equal(
-            (f1 * 1.5).freqresp([0.1, 1.0, 10])[1],
-            (h1 * 1.5).freqresp([0.1, 1.0, 10])[1])
+            (f1 * 1.5).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * 1.5).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / 1.7).freqresp([0.1, 1.0, 10])[1],
-            (h1 / 1.7).freqresp([0.1, 1.0, 10])[1])
+            (f1 / 1.7).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / 1.7).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (2.2 * f2).freqresp([0.1, 1.0, 10])[1],
-            (2.2 * h2).freqresp([0.1, 1.0, 10])[1])
+            (2.2 * f2).frequency_response([0.1, 1.0, 10])[1],
+            (2.2 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (1.3 / f2).freqresp([0.1, 1.0, 10])[1],
-            (1.3 / h2).freqresp([0.1, 1.0, 10])[1])
+            (1.3 / f2).frequency_response([0.1, 1.0, 10])[1],
+            (1.3 / h2).frequency_response([0.1, 1.0, 10])[1])
 
     def testOperatorsTf(self):
         # get two SISO transfer functions
@@ -98,24 +98,24 @@ def testOperatorsTf(self):
         f2  # reference to avoid pyflakes error
 
         np.testing.assert_array_almost_equal(
-            (f1 + h2).freqresp([0.1, 1.0, 10])[0],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 + h2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 + h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 + h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 - h2).freqresp([0.1, 1.0, 10])[0],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 - h2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 - h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 - h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[1])
         # multiplication and division
         np.testing.assert_array_almost_equal(
-            (f1 * h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 * h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 * h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 / h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 / h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / h2).frequency_response([0.1, 1.0, 10])[1])
         # the reverse does not work
 
     def testbdalg(self):
@@ -127,45 +127,45 @@ def testbdalg(self):
         f2 = FRD(h2, omega)
 
         np.testing.assert_array_almost_equal(
-            (bdalg.series(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.series(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.series(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.series(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.parallel(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.parallel(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.parallel(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.parallel(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.feedback(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.feedback(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.feedback(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.feedback(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.negate(f1)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.negate(h1)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.negate(f1)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.negate(h1)).frequency_response([0.1, 1.0, 10])[0])
 
 #       append() and connect() not implemented for FRD objects
 #        np.testing.assert_array_almost_equal(
-#            (bdalg.append(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-#            (bdalg.append(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+#            (bdalg.append(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+#            (bdalg.append(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 #
 #        f3 = bdalg.append(f1, f2, f2)
 #        h3 = bdalg.append(h1, h2, h2)
 #        Q = np.mat([ [1, 2], [2, -1] ])
 #        np.testing.assert_array_almost_equal(
-#           (bdalg.connect(f3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0],
-#            (bdalg.connect(h3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0])
+#           (bdalg.connect(f3, Q, [2], [1])).frequency_response([0.1, 1.0, 10])[0],
+#            (bdalg.connect(h3, Q, [2], [1])).frequency_response([0.1, 1.0, 10])[0])
 
     def testFeedback(self):
         h1 = TransferFunction([1], [1, 2, 2])
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(h1, omega)
         np.testing.assert_array_almost_equal(
-            f1.feedback(1).freqresp([0.1, 1.0, 10])[0],
-            h1.feedback(1).freqresp([0.1, 1.0, 10])[0])
+            f1.feedback(1).frequency_response([0.1, 1.0, 10])[0],
+            h1.feedback(1).frequency_response([0.1, 1.0, 10])[0])
 
         # Make sure default argument also works
         np.testing.assert_array_almost_equal(
-            f1.feedback().freqresp([0.1, 1.0, 10])[0],
-            h1.feedback().freqresp([0.1, 1.0, 10])[0])
+            f1.feedback().frequency_response([0.1, 1.0, 10])[0],
+            h1.feedback().frequency_response([0.1, 1.0, 10])[0])
 
     def testFeedback2(self):
         h2 = StateSpace([[-1.0, 0], [0, -2.0]], [[0.4], [0.1]],
@@ -198,11 +198,11 @@ def testMIMO(self):
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            sys.freqresp([0.1, 1.0, 10])[0],
-            f1.freqresp([0.1, 1.0, 10])[0])
+            sys.frequency_response([0.1, 1.0, 10])[0],
+            f1.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            sys.freqresp([0.1, 1.0, 10])[1],
-            f1.freqresp([0.1, 1.0, 10])[1])
+            sys.frequency_response([0.1, 1.0, 10])[1],
+            f1.frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOfb(self):
@@ -214,11 +214,11 @@ def testMIMOfb(self):
         f1 = FRD(sys, omega).feedback([[0.1, 0.3], [0.0, 1.0]])
         f2 = FRD(sys.feedback([[0.1, 0.3], [0.0, 1.0]]), omega)
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[0],
-            f2.freqresp([0.1, 1.0, 10])[0])
+            f1.frequency_response([0.1, 1.0, 10])[0],
+            f2.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[1],
-            f2.freqresp([0.1, 1.0, 10])[1])
+            f1.frequency_response([0.1, 1.0, 10])[1],
+            f2.frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOfb2(self):
@@ -232,11 +232,11 @@ def testMIMOfb2(self):
         f1 = FRD(sys, omega).feedback(K)
         f2 = FRD(sys.feedback(K), omega)
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[0],
-            f2.freqresp([0.1, 1.0, 10])[0])
+            f1.frequency_response([0.1, 1.0, 10])[0],
+            f2.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[1],
-            f2.freqresp([0.1, 1.0, 10])[1])
+            f1.frequency_response([0.1, 1.0, 10])[1],
+            f2.frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOMult(self):
@@ -248,11 +248,11 @@ def testMIMOMult(self):
         f1 = FRD(sys, omega)
         f2 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[0],
-            (sys*sys).freqresp([0.1, 1.0, 10])[0])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[0],
+            (sys*sys).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[1],
-            (sys*sys).freqresp([0.1, 1.0, 10])[1])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[1],
+            (sys*sys).frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOSmooth(self):
@@ -265,14 +265,14 @@ def testMIMOSmooth(self):
         f1 = FRD(sys, omega, smooth=True)
         f2 = FRD(sys2, omega, smooth=True)
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[0],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[0])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[0],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[1],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[1])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[1],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[2],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[2])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[2],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[2])
 
     def testAgainstOctave(self):
         # with data from octave:
@@ -285,8 +285,8 @@ def testAgainstOctave(self):
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            (f1.freqresp([1.0])[0] *
-             np.exp(1j * f1.freqresp([1.0])[1])).reshape(3, 2),
+            (f1.frequency_response([1.0])[0] *
+             np.exp(1j * f1.frequency_response([1.0])[1])).reshape(3, 2),
             np.array([[0.4 - 0.2j, 0], [0, 0.1 - 0.2j], [0, 0.3 - 0.1j]]))
 
     def test_string_representation(self, capsys):
@@ -407,13 +407,9 @@ def test_eval(self):
         with pytest.raises(ValueError):
             frd_tf.eval(2)
 
-
-    def test_evalfr_deprecated(self):
-        sys_tf = ct.tf([1], [1, 2, 1])
-        frd_tf = FRD(sys_tf, np.logspace(-1, 1, 3))
-
-        with pytest.deprecated_call():
-            frd_tf.evalfr(1.)
+        # Should get an error if we use __call__ at real-valued frequency
+        with pytest.raises(ValueError):
+            frd_tf(2)
 
     def test_repr_str(self):
         # repr printing
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index 111d4296c..e6a6934e8 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -18,7 +18,6 @@
 from control.matlab import ss, tf, bode, rss
 from control.tests.conftest import slycotonly
 
-
 pytestmark = pytest.mark.usefixtures("mplcleanup")
 
 
@@ -32,7 +31,7 @@ def test_siso():
     omega = np.linspace(10e-2, 10e2, 1000)
 
     # test frequency response
-    sys.freqresp(omega)
+    sys.frequency_response(omega)
 
     # test bode plot
     bode(sys)
@@ -97,7 +96,7 @@ def test_superimpose():
 def test_doubleint():
     """Test typcast bug with double int
 
-    30 May 2016, RMM: added to replicate typecast bug in freqresp.py
+    30 May 2016, RMM: added to replicate typecast bug in frequency_response.py
     """
     A = np.array([[0, 1], [0, 0]])
     B = np.array([[0], [1]])
@@ -117,7 +116,7 @@ def test_mimo():
     omega = np.linspace(10e-2, 10e2, 1000)
     sysMIMO = ss(A, B, C, D)
 
-    sysMIMO.freqresp(omega)
+    sysMIMO.frequency_response(omega)
     tf(sysMIMO)
 
 
@@ -215,13 +214,13 @@ def test_discrete(dsystem_type):
     omega_ok = np.linspace(10e-4, 0.99, 100) * np.pi / dsys.dt
 
     # Test frequency response
-    dsys.freqresp(omega_ok)
+    dsys.frequency_response(omega_ok)
 
     # Check for warning if frequency is out of range
     with pytest.warns(UserWarning, match="above.*Nyquist"):
         # Look for a warning about sampling above Nyquist frequency
         omega_bad = np.linspace(10e-4, 1.1, 10) * np.pi / dsys.dt
-        dsys.freqresp(omega_bad)
+        dsys.frequency_response(omega_bad)
 
     # Test bode plots (currently only implemented for SISO)
     if (dsys.inputs == 1 and dsys.outputs == 1):
@@ -332,6 +331,7 @@ def test_initial_phase(TF, initial_phase, default_phase, expected_phase):
      pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
                   -270, -3*math.pi/2, math.pi/2,    id="order5, -270"),
     ])
+
 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)
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 4a8e09930..0dcbf3325 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -1132,8 +1132,8 @@ def test_lineariosys_statespace(self, tsys):
         np.testing.assert_array_equal(
             iosys_siso.pole(), tsys.siso_linsys.pole())
         omega = np.logspace(.1, 10, 100)
-        mag_io, phase_io, omega_io = iosys_siso.freqresp(omega)
-        mag_ss, phase_ss, omega_ss = tsys.siso_linsys.freqresp(omega)
+        mag_io, phase_io, omega_io = iosys_siso.frequency_response(omega)
+        mag_ss, phase_ss, omega_ss = tsys.siso_linsys.frequency_response(omega)
         np.testing.assert_array_equal(mag_io, mag_ss)
         np.testing.assert_array_equal(phase_io, phase_ss)
         np.testing.assert_array_equal(omega_io, omega_ss)
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index 80916da1b..81a4b68b4 100644
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -93,7 +93,7 @@ def test_stability_margins_3input(tsys):
     sys, refout, refoutall = tsys
     """Test stability_margins() function with mag, phase, omega input"""
     omega = np.logspace(-2, 2, 2000)
-    mag, phase, omega_ = sys.freqresp(omega)
+    mag, phase, omega_ = sys.frequency_response(omega)
     out = stability_margins((mag, phase*180/np.pi, omega_))
     assert_allclose(out, refout, atol=1.5e-3)
 
@@ -109,7 +109,7 @@ def test_margin_3input(tsys):
     sys, refout, refoutall = tsys
     """Test margin() function with mag, phase, omega input"""
     omega = np.logspace(-2, 2, 2000)
-    mag, phase, omega_ = sys.freqresp(omega)
+    mag, phase, omega_ = sys.frequency_response(omega)
     out = margin((mag, phase*180/np.pi, omega_))
     assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
 
@@ -136,8 +136,10 @@ def test_phase_crossover_frequencies_mimo():
                            [[3], [4]]],
                           [[[1, 2, 3, 4], [1, 1]],
                            [[1, 1], [1, 1]]])
-    with pytest.raises(ControlMIMONotImplemented):
-        omega, gain = phase_crossover_frequencies(tf)
+
+    omega, gain = phase_crossover_frequencies(tf[0,0])
+    assert_allclose(omega, [1.73205,  0.], atol=1.5e-3)
+    assert_allclose(gain, [-0.5,  0.25], atol=1.5e-3)
 
 
 def test_mag_phase_omega():
@@ -145,7 +147,7 @@ def test_mag_phase_omega():
     sys = TransferFunction(15, [1, 6, 11, 6])
     out = stability_margins(sys)
     omega = np.logspace(-2, 2, 1000)
-    mag, phase, omega = sys.freqresp(omega)
+    mag, phase, omega = sys.frequency_response(omega)
     out2 = stability_margins((mag, phase*180/np.pi, omega))
     ind = [0, 1, 3, 4]   # indices of gm, pm, wg, wp -- ignore sm
     marg1 = np.array(out)[ind]
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 02fac1776..7bf90f4c6 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -341,18 +341,14 @@ def test_evalfr(self, dt, omega, resp):
         else:
             s = omega * 1j
 
-        # Correct version of the call
+        # Correct versions of the call
         np.testing.assert_allclose(evalfr(sys, s), resp, atol=1e-3)
-        # Deprecated version of the call (should generate warning)
-        with pytest.deprecated_call():
+        np.testing.assert_allclose(sys(s), resp, atol=1e-3)
+
+        # Deprecated version of the call (should generate error)
+        with pytest.raises(AttributeError):
             np.testing.assert_allclose(sys.evalfr(omega), resp, atol=1e-3)
 
-        # call above nyquist frequency
-        if dt:
-            with pytest.warns(UserWarning):
-                np.testing.assert_allclose(sys._evalfr(omega + 2 * np.pi / dt),
-                                           resp,
-                                           atol=1e-3)
 
     @slycotonly
     def test_freq_resp(self):
@@ -374,7 +370,7 @@ def test_freq_resp(self):
                       [-0.438157380501337, -1.40720969147217]]]
         true_omega = [0.1, 10.]
 
-        mag, phase, omega = sys.freqresp(true_omega)
+        mag, phase, omega = sys.frequency_response(true_omega)
 
         np.testing.assert_almost_equal(mag, true_mag)
         np.testing.assert_almost_equal(phase, true_phase)
@@ -739,9 +735,9 @@ def test_horner(self, sys322):
         sys322.horner(1. + 1.j)
 
         # Make sure result agrees with frequency response
-        mag, phase, omega = sys322.freqresp([1])
+        mag, phase, omega = sys322.frequency_response([1])
         np.testing.assert_array_almost_equal(
-            sys322.horner(1.j),
+            np.squeeze(sys322.horner(1.j)),
             mag[:, :, 0] * np.exp(1.j * phase[:, :, 0]))
 
 class TestRss:
@@ -902,7 +898,6 @@ def test_statespace_defaults(self, matarrayout):
 refkey_n = {None: 'p3', '.3g': 'p3', '.5g': 'p5'}
 refkey_r = {None: 'p', 'partitioned': 'p', 'separate': 's'}
 
-
 @pytest.mark.parametrize(" gmats,  ref",
                          [(LTX_G1, LTX_G1_REF),
                           (LTX_G2, LTX_G2_REF)])
@@ -928,3 +923,27 @@ def test_latex_repr(gmats, ref, repr_type, num_format, editsdefaults):
     refkey = "{}_{}".format(refkey_n[num_format], refkey_r[repr_type])
     assert g._repr_latex_() == ref[refkey]
 
+    def test_pole_static(self):
+        """Regression: pole() of static gain is empty array."""
+        np.testing.assert_array_equal(np.array([]),
+                                      StateSpace([], [], [], [[1]]).pole())
+
+    def test_copy_constructor(self):
+        # Create a set of matrices for a simple linear system
+        A = np.array([[-1]])
+        B = np.array([[1]])
+        C = np.array([[1]])
+        D = np.array([[0]])
+
+        # Create the first linear system and a copy
+        linsys = StateSpace(A, B, C, D)
+        cpysys = StateSpace(linsys)
+
+        # Change the original A matrix
+        A[0, 0] = -2
+        np.testing.assert_array_equal(linsys.A, [[-1]]) # original value
+        np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value
+
+        # Change the A matrix for the original system
+        linsys.A[0, 0] = -3
+        np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index c0b5e227f..1ec596a6e 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -404,34 +404,6 @@ def test_slice(self):
         assert (sys1.inputs, sys1.outputs) == (2, 1)
         assert sys1.dt == 0.5
 
-    @pytest.mark.parametrize("omega, resp",
-                             [(1, np.array([[-0.5 - 0.5j]])),
-                              (32, np.array([[0.002819593 - 0.03062847j]]))])
-    @pytest.mark.parametrize("dt", [None, 0, 1e-3])
-    def test_evalfr_siso(self, dt, omega, resp):
-        """Evaluate the frequency response at single frequencies"""
-        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
-
-        if dt:
-            sys = sample_system(sys, dt)
-            s = np.exp(omega * 1j * dt)
-        else:
-            s = omega * 1j
-
-        # Correct versions of the call
-        np.testing.assert_allclose(evalfr(sys, s), resp, atol=1e-3)
-        np.testing.assert_allclose(sys(s), resp, atol=1e-3)
-        # Deprecated version of the call (should generate warning)
-        with pytest.deprecated_call():
-            np.testing.assert_allclose(sys.evalfr(omega), resp, atol=1e-3)
-
-        # call above nyquist frequency
-        if dt:
-            with pytest.warns(UserWarning):
-                np.testing.assert_allclose(sys._evalfr(omega + 2 * np.pi / dt),
-                                           resp,
-                                           atol=1e-3)
-
     def test_is_static_gain(self):
         numstatic = 1.1
         denstatic = 1.2
@@ -454,10 +426,37 @@ def test_is_static_gain(self):
         assert not TransferFunction(numdynamicmimo,
                                     denstaticmimo).is_static_gain()
 
+    @pytest.mark.parametrize("omega, resp",
+                             [(1, np.array([[-0.5 - 0.5j]])),
+                              (32, np.array([[0.002819593 - 0.03062847j]]))])
+    @pytest.mark.parametrize("dt", [None, 0, 1e-3])
+    def test_call_siso(self, dt, omega, resp):
+        """Evaluate the frequency response of a SISO system at one frequency."""
+        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
+
+        if dt:
+            sys = sample_system(sys, dt)
+            s = np.exp(omega * 1j * dt)
+        else:
+            s = omega * 1j
+
+        # Correct versions of the call
+        np.testing.assert_allclose(evalfr(sys, s), resp, atol=1e-3)
+        np.testing.assert_allclose(sys(s), resp, atol=1e-3)
+        # Deprecated version of the call (should generate exception)
+        with pytest.raises(AttributeError):
+            np.testing.assert_allclose(sys.evalfr(omega), resp, atol=1e-3)
+
+
+    @nopython2
+    def test_call_dtime(self):
+        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1], 0.1)
+        np.testing.assert_array_almost_equal(sys(1j), -0.5 - 0.5j)
 
     @slycotonly
-    def test_evalfr_mimo(self):
-        """Evaluate the frequency response of a MIMO system at a freq"""
+    def test_call_mimo(self):
+        """Evaluate the frequency response of a MIMO system at one frequency."""
+
         num = [[[1., 2.], [0., 3.], [2., -1.]],
                [[1.], [4., 0.], [1., -4., 3.]]]
         den = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]],
@@ -467,13 +466,28 @@ def test_evalfr_mimo(self):
                 [-0.083333333333333, -0.188235294117647 - 0.847058823529412j,
                  -1. - 8.j]]
 
-        np.testing.assert_array_almost_equal(sys._evalfr(2.), resp)
+        np.testing.assert_array_almost_equal(evalfr(sys, 2j), resp)
 
         # Test call version as well
         np.testing.assert_array_almost_equal(sys(2.j), resp)
 
-    def test_freqresp_siso(self):
-        """Evaluate the SISO magnitude and phase at multiple frequencies"""
+    def test_freqresp_deprecated(self):
+        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1.])
+        # Deprecated version of the call (should generate warning)
+        import warnings
+        with warnings.catch_warnings(record=True) as w:
+            # Set up warnings filter to only show warnings in control module
+            warnings.filterwarnings("ignore")
+            warnings.filterwarnings("always", module="control")
+
+            # Make sure that we get a pending deprecation warning
+            sys.freqresp(1.)
+            assert issubclass(w[-1].category, DeprecationWarning)
+
+    def test_frequency_response_siso(self):
+        """Evaluate the magnitude and phase of a SISO system at
+        multiple frequencies."""
+
         sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
 
         truemag = [[[4.63507337473906, 0.707106781186548, 0.0866592803995351]]]
@@ -481,7 +495,7 @@ def test_freqresp_siso(self):
                        -1.32655885133871]]]
         trueomega = [0.1, 1., 10.]
 
-        mag, phase, omega = sys.freqresp(trueomega)
+        mag, phase, omega = sys.frequency_response(trueomega, squeeze=False)
 
         np.testing.assert_array_almost_equal(mag, truemag)
         np.testing.assert_array_almost_equal(phase, truephase)
@@ -509,7 +523,7 @@ def test_freqresp_mimo(self):
                        [-1.66852323, -1.89254688, -1.62050658],
                        [-0.13298964, -1.10714871, -2.75046720]]]
 
-        mag, phase, omega = sys.freqresp(true_omega)
+        mag, phase, omega = sys.frequency_response(true_omega)
 
         np.testing.assert_array_almost_equal(mag, true_mag)
         np.testing.assert_array_almost_equal(phase, true_phase)
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 9a70e36b6..a0941b74a 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -231,19 +231,70 @@ def __init__(self, *args, **kwargs):
                     dt = config.defaults['control.default_dt']
         self.dt = dt
 
-    def __call__(self, s):
-        """Evaluate the system's transfer function for a complex variable
+    def __call__(self, x, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies.
 
-        For a SISO transfer function, returns the value of the
-        transfer function.  For a MIMO transfer fuction, returns a
-        matrix of values evaluated at complex variable s."""
+        Returns the complex frequency response `sys(x)` where `x` is `s` for
+        continuous-time systems and `z` for discrete-time systems.
 
-        if self.issiso():
-            # return a scalar
-            return self.horner(s)[0][0]
+        To evaluate at a frequency omega in radians per second, enter
+        ``x = omega * 1j``, for continuous-time systems, or
+        ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use
+        :meth:`TransferFunction.frequency_response`.
+
+        Parameters
+        ----------
+        x: complex array_like or complex
+            Complex frequencies
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of `x`.
+
+        Returns
+        -------
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is `len(x)` if and
+            only if system is SISO and ``squeeze=True``.
+
+        """
+        out = self.horner(x)
+        if not hasattr(x, '__len__'):
+            # received a scalar x, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            # return a scalar/1d array of outputs
+            return out[0][0]
         else:
-            # return a matrix
-            return self.horner(s)
+            return out
+
+    def horner(self, x):
+        """Evaluate system's transfer function at complex frequency
+        using Horner's method.
+
+        Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z`
+        for discrete-time systems.
+
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
+        for a more user-friendly interface.
+
+        Parameters
+        ----------
+        x : complex array_like or complex
+            Complex frequencies
+
+        Returns
+        -------
+        output : (self.outputs, self.inputs, len(x)) complex ndarray
+            Frequency response
+
+        """
+        x_arr = np.atleast_1d(x) # force to be an array
+        out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
+        for i in range(self.outputs):
+            for j in range(self.inputs):
+                out[i][j] = (polyval(self.num[i][j], x) /
+                             polyval(self.den[i][j], x))
+        return out
 
     def _truncatecoeff(self):
         """Remove extraneous zero coefficients from num and den.
@@ -607,103 +658,19 @@ def __getitem__(self, key):
         else:
             return TransferFunction(num, den, self.dt)
 
-    def evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self._evalfr(omega) returns the value of the transfer function
-        matrix with input value s = i * omega.
-
-        """
-        warn("TransferFunction.evalfr(omega) will be deprecated in a "
-             "future release of python-control; use evalfr(sys, omega*1j) "
-             "instead", PendingDeprecationWarning)
-        return self._evalfr(omega)
-
-    def _evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency."""
-        # TODO: implement for discrete time systems
-        if isdtime(self, strict=True):
-            # Convert the frequency to discrete time
-            s = exp(1.j * omega * self.dt)
-            if np.any(omega * self.dt > pi):
-                warn("_evalfr: frequency evaluation above Nyquist frequency")
-        else:
-            s = 1.j * omega
-
-        return self.horner(s)
-
-    def horner(self, s):
-        """Evaluate the systems's transfer function for a complex variable
-
-        Returns a matrix of values evaluated at complex variable s.
-        """
-
-        # Preallocate the output.
-        if getattr(s, '__iter__', False):
-            out = empty((self.outputs, self.inputs, len(s)), dtype=complex)
-        else:
-            out = empty((self.outputs, self.inputs), dtype=complex)
-
-        for i in range(self.outputs):
-            for j in range(self.inputs):
-                out[i][j] = (polyval(self.num[i][j], s) /
-                             polyval(self.den[i][j], s))
-
-        return out
-
     def freqresp(self, omega):
-        """Evaluate the transfer function at a list of angular frequencies.
-
-        Reports the frequency response of the system,
+        """(deprecated) Evaluate transfer function at complex frequencies.
 
-             G(j*omega) = mag*exp(j*phase)
-
-        for continuous time. For discrete time systems, the response is
-        evaluated around the unit circle such that
-
-             G(exp(j*omega*dt)) = mag*exp(j*phase).
-
-        Parameters
-        ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
-
-        Returns
-        -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray or list or tuple
-            The list of sorted frequencies at which the response was
-            evaluated.
+        .. deprecated::0.9.0
+            Method has been given the more pythonic name
+            :meth:`TransferFunction.frequency_response`. Or use
+            :func:`freqresp` in the MATLAB compatibility module.
         """
-        # Preallocate outputs.
-        numfreq = len(omega)
-        mag = empty((self.outputs, self.inputs, numfreq))
-        phase = empty((self.outputs, self.inputs, numfreq))
-
-        # Figure out the frequencies
-        omega.sort()
-        if isdtime(self, strict=True):
-            slist = np.array([exp(1.j * w * self.dt) for w in omega])
-            if max(omega) * self.dt > pi:
-                warn("freqresp: frequency evaluation above Nyquist frequency")
-        else:
-            slist = np.array([1j * w for w in omega])
-
-        # Compute frequency response for each input/output pair
-        for i in range(self.outputs):
-            for j in range(self.inputs):
-                fresp = (polyval(self.num[i][j], slist) /
-                         polyval(self.den[i][j], slist))
-                mag[i, j, :] = abs(fresp)
-                phase[i, j, :] = angle(fresp)
-
-        return mag, phase, omega
+        warn("TransferFunction.freqresp(omega) will be removed in a "
+             "future release of python-control; use "
+             "sys.frequency_response(omega), or freqresp(sys, omega) in the "
+             "MATLAB compatibility module instead", DeprecationWarning)
+        return self.frequency_response(omega)
 
     def pole(self):
         """Compute the poles of a transfer function."""
diff --git a/examples/robust_mimo.py b/examples/robust_mimo.py
index d4e1335e6..9cc5a1c3b 100644
--- a/examples/robust_mimo.py
+++ b/examples/robust_mimo.py
@@ -43,7 +43,7 @@ def triv_sigma(g, w):
     g - LTI object, order n
     w - frequencies, length m
     s - (m,n) array of singular values of g(1j*w)"""
-    m, p, _ = g.freqresp(w)
+    m, p, _ = g.frequency_response(w)
     sjw = (m*np.exp(1j*p)).transpose(2, 0, 1)
     sv = np.linalg.svd(sjw, compute_uv=False)
     return sv
diff --git a/examples/robust_siso.py b/examples/robust_siso.py
index 87fcdb707..17ce10927 100644
--- a/examples/robust_siso.py
+++ b/examples/robust_siso.py
@@ -50,10 +50,10 @@
 
 # frequency response
 omega = np.logspace(-2, 2, 101)
-ws1mag, _, _ = ws1.freqresp(omega)
-s1mag, _, _ = s1.freqresp(omega)
-ws2mag, _, _ = ws2.freqresp(omega)
-s2mag, _, _ = s2.freqresp(omega)
+ws1mag, _, _ = ws1.frequency_response(omega)
+s1mag, _, _ = s1.frequency_response(omega)
+ws2mag, _, _ = ws2.frequency_response(omega)
+s2mag, _, _ = s2.frequency_response(omega)
 
 plt.figure(1)
 # text uses log-scaled absolute, but dB are probably more familiar to most control engineers