diff --git a/control/iosys.py b/control/iosys.py
index 479039c3d..1b55053e3 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -32,7 +32,8 @@
 from warnings import warn
 
 from .statesp import StateSpace, tf2ss, _convert_to_statespace
-from .timeresp import _check_convert_array, _process_time_response
+from .timeresp import _check_convert_array, _process_time_response, \
+    TimeResponseData
 from .lti import isctime, isdtime, common_timebase
 from . import config
 
@@ -878,7 +879,7 @@ def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None,
     def __call__(sys, u, params=None, squeeze=None):
         """Evaluate a (static) nonlinearity at a given input value
 
-        If a nonlinear I/O system has not internal state, then evaluating the
+        If a nonlinear I/O system has no internal state, then evaluating the
         system at an input `u` gives the output `y = F(u)`, determined by the
         output function.
 
@@ -907,7 +908,8 @@ def __call__(sys, u, params=None, squeeze=None):
 
         # Evaluate the function on the argument
         out = sys._out(0, np.array((0,)), np.asarray(u))
-        _, out = _process_time_response(sys, None, out, None, squeeze=squeeze)
+        _, out = _process_time_response(
+            None, out, issiso=sys.issiso(), squeeze=squeeze)
         return out
 
     def _update_params(self, params, warning=False):
@@ -1484,14 +1486,21 @@ def input_output_response(
     ----------
     sys : InputOutputSystem
         Input/output system to simulate.
+
     T : array-like
         Time steps at which the input is defined; values must be evenly spaced.
+
     U : array-like or number, optional
         Input array giving input at each time `T` (default = 0).
+
     X0 : array-like or number, optional
         Initial condition (default = 0).
+
     return_x : bool, optional
+        If True, return the state vector when assigning to a tuple (default =
+        False).  See :func:`forced_response` for more details.
         If True, return the values of the state at each time (default = False).
+
     squeeze : bool, optional
         If True and if the system has a single output, return the system
         output as a 1D array rather than a 2D array.  If False, return the
@@ -1500,15 +1509,27 @@ def input_output_response(
 
     Returns
     -------
-    T : array
-        Time values of the output.
-    yout : array
-        Response of the system.  If the system is SISO and squeeze is not
-        True, the array is 1D (indexed by time).  If the system is not SISO or
-        squeeze is False, the array is 2D (indexed by the output number and
-        time).
-    xout : array
-        Time evolution of the state vector (if return_x=True).
+    results : TimeResponseData
+        Time response represented as a :class:`TimeResponseData` object
+        containing the following properties:
+
+        * time (array): Time values of the output.
+
+        * outputs (array): Response of the system.  If the system is SISO and
+          `squeeze` is not True, the array is 1D (indexed by time).  If the
+          system is not SISO or `squeeze` is False, the array is 2D (indexed
+          by output and time).
+
+        * states (array): Time evolution of the state vector, represented as
+          a 2D array indexed by state and time.
+
+        * inputs (array): Input(s) to the system, indexed by input and time.
+
+        The return value of the system can also be accessed by assigning the
+        function to a tuple of length 2 (time, output) or of length 3 (time,
+        output, state) if ``return_x`` is ``True``.  If the input/output
+        system signals are named, these names will be used as labels for the
+        time response.
 
     Other parameters
     ----------------
@@ -1570,8 +1591,9 @@ def input_output_response(
         for i in range(len(T)):
             u = U[i] if len(U.shape) == 1 else U[:, i]
             y[:, i] = sys._out(T[i], [], u)
-        return _process_time_response(
-            sys, T, y, np.array((0, 0, np.asarray(T).size)),
+        return TimeResponseData(
+            T, y, None, U, issiso=sys.issiso(),
+            output_labels=sys.output_index, input_labels=sys.input_index,
             transpose=transpose, return_x=return_x, squeeze=squeeze)
 
     # create X0 if not given, test if X0 has correct shape
@@ -1666,8 +1688,11 @@ def ivp_rhs(t, x):
     else:                       # Neither ctime or dtime??
         raise TypeError("Can't determine system type")
 
-    return _process_time_response(sys, soln.t, y, soln.y, transpose=transpose,
-                                  return_x=return_x, squeeze=squeeze)
+    return TimeResponseData(
+        soln.t, y, soln.y, U, issiso=sys.issiso(),
+        output_labels=sys.output_index, input_labels=sys.input_index,
+        state_labels=sys.state_index,
+        transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
 def find_eqpt(sys, x0, u0=[], y0=None, t=0, params={},
diff --git a/control/optimal.py b/control/optimal.py
index b88513f69..76e9a2d31 100644
--- a/control/optimal.py
+++ b/control/optimal.py
@@ -16,7 +16,7 @@
 import logging
 import time
 
-from .timeresp import _process_time_response
+from .timeresp import TimeResponseData
 
 __all__ = ['find_optimal_input']
 
@@ -826,13 +826,14 @@ def __init__(
         else:
             states = None
 
-        retval = _process_time_response(
-            ocp.system, ocp.timepts, inputs, states,
+        # Process data as a time response (with "outputs" = inputs)
+        response = TimeResponseData(
+            ocp.timepts, inputs, states, issiso=ocp.system.issiso(),
             transpose=transpose, return_x=return_states, squeeze=squeeze)
 
-        self.time = retval[0]
-        self.inputs = retval[1]
-        self.states = None if states is None else retval[2]
+        self.time = response.time
+        self.inputs = response.outputs
+        self.states = response.states
 
 
 # Compute the input for a nonlinear, (constrained) optimal control problem
diff --git a/control/statesp.py b/control/statesp.py
index 6b3a1dff3..6a46a26e0 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -1932,10 +1932,10 @@ def rss(states=1, outputs=1, inputs=1, strictly_proper=False):
     ----------
     states : int
         Number of state variables
-    inputs : int
-        Number of system inputs
     outputs : int
         Number of system outputs
+    inputs : int
+        Number of system inputs
     strictly_proper : bool, optional
         If set to 'True', returns a proper system (no direct term).
 
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index c74c0c06d..435d8a60c 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -1117,7 +1117,7 @@ def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape1, shape2):
     @pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.ss2io])
     def test_squeeze_exception(self, fcn):
         sys = fcn(ct.rss(2, 1, 1))
-        with pytest.raises(ValueError, match="unknown squeeze value"):
+        with pytest.raises(ValueError, match="Unknown squeeze value"):
             step_response(sys, squeeze=1)
 
     @pytest.mark.usefixtures("editsdefaults")
diff --git a/control/tests/trdata_test.py b/control/tests/trdata_test.py
new file mode 100644
index 000000000..fcd8676e9
--- /dev/null
+++ b/control/tests/trdata_test.py
@@ -0,0 +1,362 @@
+"""trdata_test.py - test return values from time response functions
+
+RMM, 22 Aug 2021
+
+This set of unit tests covers checks to make sure that the various time
+response functions are returning the right sets of objects in the (new)
+InputOutputResponse class.
+
+"""
+
+import pytest
+
+import numpy as np
+import control as ct
+
+
+@pytest.mark.parametrize(
+    "nout, nin, squeeze", [
+        [1,   1,  None],
+        [1,   1,  True],
+        [1,   1,  False],
+        [1,   2,  None],
+        [1,   2,  True],
+        [1,   2,  False],
+        [2,   1,  None],
+        [2,   1,  True],
+        [2,   1,  False],
+        [2,   3,  None],
+        [2,   3,  True],
+        [2,   3,  False],
+])
+def test_trdata_shapes(nin, nout, squeeze):
+    # SISO, single trace
+    sys = ct.rss(4, nout, nin, strictly_proper=True)
+    T = np.linspace(0, 1, 10)
+    U = np.outer(np.ones(nin), np.sin(T) )
+    X0 = np.ones(sys.nstates)
+
+    #
+    # Initial response
+    #
+    res = ct.initial_response(sys, X0=X0)
+    ntimes = res.time.shape[0]
+
+    # Check shape of class members
+    assert len(res.time.shape) == 1
+    assert res.y.shape == (sys.noutputs, ntimes)
+    assert res.x.shape == (sys.nstates, ntimes)
+    assert res.u is None
+
+    # Check dimensions of the response
+    assert res.ntraces == 0     # single trace
+    assert res.ninputs == 0     # no input for initial response
+    assert res.noutputs == sys.noutputs
+    assert res.nstates == sys.nstates
+
+    # Check shape of class properties
+    if sys.issiso():
+        assert res.outputs.shape == (ntimes,)
+        assert res._legacy_states.shape == (sys.nstates, ntimes)
+        assert res.states.shape == (sys.nstates, ntimes)
+        assert res.inputs is None
+    elif res.squeeze is True:
+        assert res.outputs.shape == (ntimes, )
+        assert res._legacy_states.shape == (sys.nstates, ntimes)
+        assert res.states.shape == (sys.nstates, ntimes)
+        assert res.inputs is None
+    else:
+        assert res.outputs.shape == (sys.noutputs, ntimes)
+        assert res._legacy_states.shape == (sys.nstates, ntimes)
+        assert res.states.shape == (sys.nstates, ntimes)
+        assert res.inputs is None
+
+    #
+    # Impulse and step response
+    #
+    for fcn in (ct.impulse_response, ct.step_response):
+        res = fcn(sys, squeeze=squeeze)
+        ntimes = res.time.shape[0]
+
+        # Check shape of class members
+        assert len(res.time.shape) == 1
+        assert res.y.shape == (sys.noutputs, sys.ninputs, ntimes)
+        assert res.x.shape == (sys.nstates, sys.ninputs, ntimes)
+        assert res.u.shape == (sys.ninputs, sys.ninputs, ntimes)
+
+        # Check shape of class members
+        assert res.ntraces == sys.ninputs
+        assert res.ninputs == sys.ninputs
+        assert res.noutputs == sys.noutputs
+        assert res.nstates == sys.nstates
+
+        # Check shape of inputs and outputs
+        if sys.issiso() and squeeze is not False:
+            assert res.outputs.shape == (ntimes, )
+            assert res.states.shape == (sys.nstates, ntimes)
+            assert res.inputs.shape == (ntimes, )
+        elif res.squeeze is True:
+            assert res.outputs.shape == \
+                np.empty((sys.noutputs, sys.ninputs, ntimes)).squeeze().shape
+            assert res.states.shape == \
+                np.empty((sys.nstates, sys.ninputs, ntimes)).squeeze().shape
+            assert res.inputs.shape == \
+                np.empty((sys.ninputs, sys.ninputs, ntimes)).squeeze().shape
+        else:
+            assert res.outputs.shape == (sys.noutputs, sys.ninputs, ntimes)
+            assert res.states.shape == (sys.nstates, sys.ninputs, ntimes)
+            assert res.inputs.shape == (sys.ninputs, sys.ninputs, ntimes)
+
+        # Check legacy state space dimensions (not affected by squeeze)
+        if sys.issiso():
+            assert res._legacy_states.shape == (sys.nstates, ntimes)
+        else:
+            assert res._legacy_states.shape == \
+                (sys.nstates, sys.ninputs, ntimes)
+
+    #
+    # Forced response
+    #
+    res = ct.forced_response(sys, T, U, X0, squeeze=squeeze)
+    ntimes = res.time.shape[0]
+
+    # Check shape of class members
+    assert len(res.time.shape) == 1
+    assert res.y.shape == (sys.noutputs, ntimes)
+    assert res.x.shape == (sys.nstates, ntimes)
+    assert res.u.shape == (sys.ninputs, ntimes)
+
+    # Check dimensions of the response
+    assert res.ntraces == 0     # single trace
+    assert res.ninputs == sys.ninputs
+    assert res.noutputs == sys.noutputs
+    assert res.nstates == sys.nstates
+
+    # Check shape of inputs and outputs
+    if sys.issiso() and squeeze is not False:
+        assert res.outputs.shape == (ntimes,)
+        assert res.states.shape == (sys.nstates, ntimes)
+        assert res.inputs.shape == (ntimes,)
+    elif squeeze is True:
+        assert res.outputs.shape == \
+            np.empty((sys.noutputs, 1, ntimes)).squeeze().shape
+        assert res.states.shape == \
+            np.empty((sys.nstates, 1, ntimes)).squeeze().shape
+        assert res.inputs.shape == \
+            np.empty((sys.ninputs, 1, ntimes)).squeeze().shape
+    else:                       # MIMO or squeeze is False
+        assert res.outputs.shape == (sys.noutputs, ntimes)
+        assert res.states.shape == (sys.nstates, ntimes)
+        assert res.inputs.shape == (sys.ninputs, ntimes)
+
+    # Check state space dimensions (not affected by squeeze)
+    assert res.states.shape == (sys.nstates, ntimes)
+
+
+def test_response_copy():
+    # Generate some initial data to use
+    sys_siso = ct.rss(4, 1, 1)
+    response_siso = ct.step_response(sys_siso)
+    siso_ntimes = response_siso.time.size
+
+    sys_mimo = ct.rss(4, 2, 1)
+    response_mimo = ct.step_response(sys_mimo)
+    mimo_ntimes = response_mimo.time.size
+
+    # Transpose
+    response_mimo_transpose = response_mimo(transpose=True)
+    assert response_mimo.outputs.shape == (2, 1, mimo_ntimes)
+    assert response_mimo_transpose.outputs.shape == (mimo_ntimes, 2, 1)
+    assert response_mimo.states.shape == (4, 1, mimo_ntimes)
+    assert response_mimo_transpose.states.shape == (mimo_ntimes, 4, 1)
+
+    # Squeeze
+    response_siso_as_mimo = response_siso(squeeze=False)
+    assert response_siso_as_mimo.outputs.shape == (1, 1, siso_ntimes)
+    assert response_siso_as_mimo.states.shape == (4, 1, siso_ntimes)
+    assert response_siso_as_mimo._legacy_states.shape == (4, siso_ntimes)
+
+    response_mimo_squeezed = response_mimo(squeeze=True)
+    assert response_mimo_squeezed.outputs.shape == (2, mimo_ntimes)
+    assert response_mimo_squeezed.states.shape == (4, mimo_ntimes)
+    assert response_mimo_squeezed._legacy_states.shape == (4, 1, mimo_ntimes)
+
+    # Squeeze and transpose
+    response_mimo_sqtr = response_mimo(squeeze=True, transpose=True)
+    assert response_mimo_sqtr.outputs.shape == (mimo_ntimes, 2)
+    assert response_mimo_sqtr.states.shape == (mimo_ntimes, 4)
+    assert response_mimo_sqtr._legacy_states.shape == (mimo_ntimes, 4, 1)
+
+    # Return_x
+    t, y = response_mimo
+    t, y = response_mimo()
+    t, y, x = response_mimo(return_x=True)
+    with pytest.raises(ValueError, match="too many"):
+        t, y = response_mimo(return_x=True)
+    with pytest.raises(ValueError, match="not enough"):
+        t, y, x = response_mimo
+
+    # Labels
+    assert response_mimo.output_labels is None
+    assert response_mimo.state_labels is None
+    assert response_mimo.input_labels is None
+    response = response_mimo(
+        output_labels=['y1', 'y2'], input_labels='u',
+        state_labels=["x[%d]" % i for i in range(4)])
+    assert response.output_labels == ['y1', 'y2']
+    assert response.state_labels == ['x[0]', 'x[1]', 'x[2]', 'x[3]']
+    assert response.input_labels == ['u']
+
+    # Unknown keyword
+    with pytest.raises(ValueError, match="Unknown parameter(s)*"):
+        response_bad_kw = response_mimo(input=0)
+
+
+def test_trdata_labels():
+    # Create an I/O system with labels
+    sys = ct.rss(4, 3, 2)
+    iosys = ct.LinearIOSystem(sys)
+
+    T = np.linspace(1, 10, 10)
+    U = [np.sin(T), np.cos(T)]
+
+    # Create a response
+    response = ct.input_output_response(iosys, T, U)
+
+    # Make sure the labels got created
+    np.testing.assert_equal(
+        response.output_labels, ["y[%d]" % i for i in range(sys.noutputs)])
+    np.testing.assert_equal(
+        response.state_labels, ["x[%d]" % i for i in range(sys.nstates)])
+    np.testing.assert_equal(
+        response.input_labels, ["u[%d]" % i for i in range(sys.ninputs)])
+
+
+def test_trdata_multitrace():
+    #
+    # Output signal processing
+    #
+
+    # Proper call of multi-trace data w/ ambiguous 2D output
+    response = ct.TimeResponseData(
+        np.zeros(5), np.ones((2, 5)), np.zeros((3, 2, 5)),
+        np.ones((4, 2, 5)), multi_trace=True)
+    assert response.ntraces == 2
+    assert response.noutputs == 1
+    assert response.nstates == 3
+    assert response.ninputs == 4
+
+    # Proper call of single trace w/ ambiguous 2D output
+    response = ct.TimeResponseData(
+        np.zeros(5), np.ones((2, 5)), np.zeros((3, 5)),
+        np.ones((4, 5)), multi_trace=False)
+    assert response.ntraces == 0
+    assert response.noutputs == 2
+    assert response.nstates == 3
+    assert response.ninputs == 4
+
+    # Proper call of multi-trace data w/ ambiguous 1D output
+    response = ct.TimeResponseData(
+        np.zeros(5), np.ones(5), np.zeros((3, 5)),
+        np.ones((4, 5)), multi_trace=False)
+    assert response.ntraces == 0
+    assert response.noutputs == 1
+    assert response.nstates == 3
+    assert response.ninputs == 4
+    assert response.y.shape == (1, 5)           # Make sure reshape occured
+
+    # Output vector not the right shape
+    with pytest.raises(ValueError, match="Output vector is the wrong shape"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones((1, 2, 3, 5)), None, None)
+
+    # Inconsistent output vector: different number of time points
+    with pytest.raises(ValueError, match="Output vector does not match time"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones(6), np.zeros(5), np.zeros(5))
+
+    #
+    # State signal processing
+    #
+
+    # For multi-trace, state must be 3D
+    with pytest.raises(ValueError, match="State vector is the wrong shape"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones((1, 5)), np.zeros((3, 5)), multi_trace=True)
+
+    # If not multi-trace, state must be 2D
+    with pytest.raises(ValueError, match="State vector is the wrong shape"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones(5), np.zeros((3, 1, 5)), multi_trace=False)
+
+    # State vector in the wrong shape
+    with pytest.raises(ValueError, match="State vector is the wrong shape"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones((1, 2, 5)), np.zeros((2, 1, 5)))
+
+    # Inconsistent state vector: different number of time points
+    with pytest.raises(ValueError, match="State vector does not match time"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones(5), np.zeros((1, 6)), np.zeros(5))
+
+    #
+    # Input signal processing
+    #
+
+    # Proper call of multi-trace data with 2D input
+    response = ct.TimeResponseData(
+        np.zeros(5), np.ones((2, 5)), np.zeros((3, 2, 5)),
+        np.ones((2, 5)), multi_trace=True)
+    assert response.ntraces == 2
+    assert response.noutputs == 1
+    assert response.nstates == 3
+    assert response.ninputs == 1
+
+    # Input vector in the wrong shape
+    with pytest.raises(ValueError, match="Input vector is the wrong shape"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones((1, 2, 5)), None, np.zeros((2, 1, 5)))
+
+    # Inconsistent input vector: different number of time points
+    with pytest.raises(ValueError, match="Input vector does not match time"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones(5), np.zeros((1, 5)), np.zeros(6))
+
+
+def test_trdata_exceptions():
+    # Incorrect dimension for time vector
+    with pytest.raises(ValueError, match="Time vector must be 1D"):
+        ct.TimeResponseData(np.zeros((2,2)), np.zeros(2), None)
+
+    # Infer SISO system from inputs and outputs
+    response = ct.TimeResponseData(
+        np.zeros(5), np.ones(5), None, np.ones(5))
+    assert response.issiso
+
+    response = ct.TimeResponseData(
+        np.zeros(5), np.ones((1, 5)), None, np.ones((1, 5)))
+    assert response.issiso
+
+    response = ct.TimeResponseData(
+        np.zeros(5), np.ones((1, 2, 5)), None, np.ones((1, 2, 5)))
+    assert response.issiso
+
+    # Not enough input to infer whether SISO
+    with pytest.raises(ValueError, match="Can't determine if system is SISO"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones((1, 2, 5)), np.ones((4, 2, 5)), None)
+
+    # Not enough input to infer whether SISO
+    with pytest.raises(ValueError, match="Keyword `issiso` does not match"):
+        response = ct.TimeResponseData(
+            np.zeros(5), np.ones((2, 5)), None, np.ones((1, 5)), issiso=True)
+
+    # Unknown squeeze keyword value
+    with pytest.raises(ValueError, match="Unknown squeeze value"):
+        response=ct.TimeResponseData(
+            np.zeros(5), np.ones(5), None, np.ones(5), squeeze=1)
+
+    # Legacy interface index error
+    response[0], response[1], response[2]
+    with pytest.raises(IndexError):
+        response[3]
diff --git a/control/timeresp.py b/control/timeresp.py
index 989a832cb..75e1dcf0b 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -64,6 +64,9 @@
 Modified by Ilhan Polat to improve automatic time vector creation
 Date: August 17, 2020
 
+Modified by Richard Murray to add TimeResponseData class
+Date: August 2021
+
 $Id$
 """
 
@@ -73,19 +76,621 @@
 import scipy as sp
 from numpy import einsum, maximum, minimum
 from scipy.linalg import eig, eigvals, matrix_balance, norm
+from copy import copy
 
 from . import config
 from .lti import isctime, isdtime
 from .statesp import StateSpace, _convert_to_statespace, _mimo2simo, _mimo2siso
 from .xferfcn import TransferFunction
 
-__all__ = ['forced_response', 'step_response', 'step_info', 'initial_response',
-           'impulse_response']
+__all__ = ['forced_response', 'step_response', 'step_info',
+           'initial_response', 'impulse_response', 'TimeResponseData']
+
+
+class TimeResponseData():
+    """A class for returning time responses.
+
+    This class maintains and manipulates the data corresponding to the
+    temporal response of an input/output system.  It is used as the return
+    type for time domain simulations (step response, input/output response,
+    etc).
+
+    A time response consists of a time vector, an output vector, and
+    optionally an input vector and/or state vector.  Inputs and outputs can
+    be 1D (scalar input/output) or 2D (vector input/output).
+
+    A time response can be stored for multiple input signals (called traces),
+    with the output and state indexed by the trace number.  This allows for
+    input/output response matrices, which is mainly useful for impulse and
+    step responses for linear systems.  For multi-trace responses, the same
+    time vector must be used for all traces.
+
+    Time responses are accessed through either the raw data, stored as
+    :attr:`t`, :attr:`y`, :attr:`x`, :attr:`u`, or using a set of properties
+    :attr:`time`, :attr:`outputs`, :attr:`states`, :attr:`inputs`.  When
+    accessing time responses via their properties, squeeze processing is
+    applied so that (by default) single-input, single-output systems will have
+    the output and input indices supressed.  This behavior is set using the
+    ``squeeze`` keyword.
+
+    Attributes
+    ----------
+    t : 1D array
+        Time values of the input/output response(s).  This attribute is
+        normally accessed via the :attr:`time` property.
+
+    y : 2D or 3D array
+        Output response data, indexed either by output index and time (for
+        single trace responses) or output, trace, and time (for multi-trace
+        responses).  These data are normally accessed via the :attr:`outputs`
+        property, which performs squeeze processing.
+
+    x : 2D or 3D array, or None
+        State space data, indexed either by output number and time (for single
+        trace responses) or output, trace, and time (for multi-trace
+        responses).  If no state data are present, value is ``None``. These
+        data are normally accessed via the :attr:`states` property, which
+        performs squeeze processing.
+
+    u : 2D or 3D array, or None
+        Input signal data, indexed either by input index and time (for single
+        trace responses) or input, trace, and time (for multi-trace
+        responses).  If no input data are present, value is ``None``.  These
+        data are normally accessed via the :attr:`inputs` property, which
+        performs squeeze processing.
+
+    squeeze : bool, optional
+        By default, if a system is single-input, single-output (SISO)
+        then the outputs (and inputs) are returned as a 1D array
+        (indexed by time) and if a system is multi-input or
+        multi-output, then the outputs are returned as a 2D array
+        (indexed by output and time) or a 3D array (indexed by output,
+        trace, and time).  If ``squeeze=True``, access to the output
+        response will remove single-dimensional entries from the shape
+        of the inputs and outputs even if the system is not SISO. If
+        ``squeeze=False``, the output is returned as a 2D or 3D array
+        (indexed by the output [if multi-input], trace [if multi-trace]
+        and time) even if the system is SISO. The default value can be
+        set using config.defaults['control.squeeze_time_response'].
+
+    transpose : bool, optional
+        If True, transpose all input and output arrays (for backward
+        compatibility with MATLAB and :func:`scipy.signal.lsim`).  Default
+        value is False.
+
+    issiso : bool, optional
+        Set to ``True`` if the system generating the data is single-input,
+        single-output.  If passed as ``None`` (default), the input data
+        will be used to set the value.
+
+    ninputs, noutputs, nstates : int
+        Number of inputs, outputs, and states of the underlying system.
+
+    input_labels, output_labels, state_labels : array of str
+        Names for the input, output, and state variables.
+
+    ntraces : int
+        Number of independent traces represented in the input/output
+        response.  If ntraces is 0 then the data represents a single trace
+        with the trace index surpressed in the data.
+
+    Notes
+    -----
+    1. For backward compatibility with earlier versions of python-control,
+       this class has an ``__iter__`` method that allows it to be assigned
+       to a tuple with a variable number of elements.  This allows the
+       following patterns to work:
+
+         t, y = step_response(sys)
+         t, y, x = step_response(sys, return_x=True)
+
+       When using this (legacy) interface, the state vector is not affected by
+       the `squeeze` parameter.
+
+    2. For backward compatibility with earlier version of python-control,
+       this class has ``__getitem__`` and ``__len__`` methods that allow the
+       return value to be indexed:
+
+         response[0]: returns the time vector
+         response[1]: returns the output vector
+         response[2]: returns the state vector
+
+       When using this (legacy) interface, the state vector is not affected by
+       the `squeeze` parameter.
+
+    3. The default settings for ``return_x``, ``squeeze`` and ``transpose``
+       can be changed by calling the class instance and passing new values:
+
+         response(tranpose=True).input
+
+       See :meth:`TimeResponseData.__call__` for more information.
+
+    """
+
+    def __init__(
+            self, time, outputs, states=None, inputs=None, issiso=None,
+            output_labels=None, state_labels=None, input_labels=None,
+            transpose=False, return_x=False, squeeze=None, multi_trace=False
+    ):
+        """Create an input/output time response object.
+
+        Parameters
+        ----------
+        time : 1D array
+            Time values of the output.  Ignored if None.
+
+        outputs : ndarray
+            Output response of the system.  This can either be a 1D array
+            indexed by time (for SISO systems or MISO systems with a specified
+            input), a 2D array indexed by output and time (for MIMO systems
+            with no input indexing, such as initial_response or forced
+            response) or trace and time (for SISO systems with multiple
+            traces), or a 3D array indexed by output, trace, and time (for
+            multi-trace input/output responses).
+
+        states : array, optional
+            Individual response of each state variable. This should be a 2D
+            array indexed by the state index and time (for single trace
+            systems) or a 3D array indexed by state, trace, and time.
+
+        inputs : array, optional
+            Inputs used to generate the output.  This can either be a 1D
+            array indexed by time (for SISO systems or MISO/MIMO systems
+            with a specified input), a 2D array indexed either by input and
+            time (for a multi-input system) or trace and time (for a
+            single-input, multi-trace response), or a 3D array indexed by
+            input, trace, and time.
+
+        sys : LTI or InputOutputSystem, optional
+            System that generated the data.  If desired, the system used to
+            generate the data can be stored along with the data.
+
+        squeeze : bool, optional
+            By default, if a system is single-input, single-output (SISO)
+            then the inputs and outputs are returned as a 1D array (indexed
+            by time) and if a system is multi-input or multi-output, then
+            the inputs are returned as a 2D array (indexed by input and
+            time) and the outputs are returned as either a 2D array (indexed
+            by output and time) or a 3D array (indexed by output, trace, and
+            time).  If squeeze=True, access to the output response will
+            remove single-dimensional entries from the shape of the inputs
+            and outputs even if the system is not SISO. If squeeze=False,
+            keep the input as a 2D or 3D array (indexed by the input (if
+            multi-input), trace (if single input) and time) and the output
+            as a 3D array (indexed by the output, trace, and time) even if
+            the system is SISO. The default value can be set using
+            config.defaults['control.squeeze_time_response'].
+
+        Other parameters
+        ----------------
+        input_labels, output_labels, state_labels: array of str, optional
+            Optional labels for the inputs, outputs, and states, given as a
+            list of strings matching the appropriate signal dimension.
+
+        transpose : bool, optional
+            If True, transpose all input and output arrays (for backward
+            compatibility with MATLAB and :func:`scipy.signal.lsim`).
+            Default value is False.
+
+        return_x : bool, optional
+            If True, return the state vector when enumerating result by
+            assigning to a tuple (default = False).
+
+        multi_trace : bool, optional
+            If ``True``, then 2D input array represents multiple traces.  For
+            a MIMO system, the ``input`` attribute should then be set to
+            indicate which trace is being specified.  Default is ``False``.
+
+        """
+        #
+        # Process and store the basic input/output elements
+        #
+
+        # Time vector
+        self.t = np.atleast_1d(time)
+        if self.t.ndim != 1:
+            raise ValueError("Time vector must be 1D array")
+
+        #
+        # Output vector (and number of traces)
+        #
+        self.y = np.array(outputs)
+
+        if self.y.ndim == 3:
+            multi_trace = True
+            self.noutputs = self.y.shape[0]
+            self.ntraces = self.y.shape[1]
+
+        elif multi_trace and self.y.ndim == 2:
+            self.noutputs = 1
+            self.ntraces = self.y.shape[0]
+
+        elif not multi_trace and self.y.ndim == 2:
+            self.noutputs = self.y.shape[0]
+            self.ntraces = 0
+
+        elif not multi_trace and self.y.ndim == 1:
+            self.noutputs = 1
+            self.ntraces = 0
+
+            # Reshape the data to be 2D for consistency
+            self.y = self.y.reshape(self.noutputs, -1)
+
+        else:
+            raise ValueError("Output vector is the wrong shape")
+
+        # Check and store labels, if present
+        self.output_labels = _process_labels(
+            output_labels, "output", self.noutputs)
+
+        # Make sure time dimension of output is the right length
+        if self.t.shape[-1] != self.y.shape[-1]:
+            raise ValueError("Output vector does not match time vector")
+
+        #
+        # State vector (optional)
+        #
+        # If present, the shape of the state vector should be consistent
+        # with the multi-trace nature of the data.
+        #
+        if states is None:
+            self.x = None
+            self.nstates = 0
+        else:
+            self.x = np.array(states)
+            self.nstates = self.x.shape[0]
+
+            # Make sure the shape is OK
+            if multi_trace and \
+               (self.x.ndim != 3 or self.x.shape[1] != self.ntraces) or \
+               not multi_trace and self.x.ndim != 2 :
+                raise ValueError("State vector is the wrong shape")
+
+            # Make sure time dimension of state is the right length
+            if self.t.shape[-1] != self.x.shape[-1]:
+                raise ValueError("State vector does not match time vector")
+
+        # Check and store labels, if present
+        self.state_labels = _process_labels(
+            state_labels, "state", self.nstates)
+
+        #
+        # Input vector (optional)
+        #
+        # If present, the shape and dimensions of the input vector should be
+        # consistent with the trace count computed above.
+        #
+        if inputs is None:
+            self.u = None
+            self.ninputs = 0
+
+        else:
+            self.u = np.array(inputs)
+
+            # Make sure the shape is OK and figure out the nuumber of inputs
+            if multi_trace and self.u.ndim == 3 and \
+               self.u.shape[1] == self.ntraces:
+                self.ninputs = self.u.shape[0]
+
+            elif multi_trace and self.u.ndim == 2 and \
+                 self.u.shape[0] == self.ntraces:
+                self.ninputs = 1
+
+            elif not multi_trace and self.u.ndim == 2 and \
+                 self.ntraces == 0:
+                self.ninputs = self.u.shape[0]
+
+            elif not multi_trace and self.u.ndim == 1:
+                self.ninputs = 1
+
+                # Reshape the data to be 2D for consistency
+                self.u = self.u.reshape(self.ninputs, -1)
+
+            else:
+                raise ValueError("Input vector is the wrong shape")
+
+            # Make sure time dimension of output is the right length
+            if self.t.shape[-1] != self.u.shape[-1]:
+                raise ValueError("Input vector does not match time vector")
+
+        # Check and store labels, if present
+        self.input_labels = _process_labels(
+            input_labels, "input", self.ninputs)
+
+        # Figure out if the system is SISO
+        if issiso is None:
+            # Figure out based on the data
+            if self.ninputs == 1:
+                issiso = (self.noutputs == 1)
+            elif self.ninputs > 1:
+                issiso = False
+            else:
+                # Missing input data => can't resolve
+                raise ValueError("Can't determine if system is SISO")
+        elif issiso is True and (self.ninputs > 1 or self.noutputs > 1):
+            raise ValueError("Keyword `issiso` does not match data")
+
+        # Set the value to be used for future processing
+        self.issiso = issiso
+
+        # Keep track of whether to squeeze inputs, outputs, and states
+        if not (squeeze is True or squeeze is None or squeeze is False):
+            raise ValueError("Unknown squeeze value")
+        self.squeeze = squeeze
+
+        # Keep track of whether to transpose for MATLAB/scipy.signal
+        self.transpose = transpose
+
+        # Store legacy keyword values (only needed for legacy interface)
+        self.return_x = return_x
+
+    def __call__(self, **kwargs):
+        """Change value of processing keywords.
+
+        Calling the time response object will create a copy of the object and
+        change the values of the keywords used to control the ``outputs``,
+        ``states``, and ``inputs`` properties.
+
+        Parameters
+        ----------
+        squeeze : bool, optional
+            If squeeze=True, access to the output response will remove
+            single-dimensional entries from the shape of the inputs, outputs,
+            and states even if the system is not SISO. If squeeze=False, keep
+            the input as a 2D or 3D array (indexed by the input (if
+            multi-input), trace (if single input) and time) and the output and
+            states as a 3D array (indexed by the output/state, trace, and
+            time) even if the system is SISO.
+
+        transpose : bool, optional
+            If True, transpose all input and output arrays (for backward
+            compatibility with MATLAB and :func:`scipy.signal.lsim`).
+            Default value is False.
+
+        return_x : bool, optional
+            If True, return the state vector when enumerating result by
+            assigning to a tuple (default = False).
+
+        input_labels, output_labels, state_labels: array of str
+            Labels for the inputs, outputs, and states, given as a
+            list of strings matching the appropriate signal dimension.
+
+        """
+        # Make a copy of the object
+        response = copy(self)
+
+        # Update any keywords that we were passed
+        response.transpose = kwargs.pop('transpose', self.transpose)
+        response.squeeze = kwargs.pop('squeeze', self.squeeze)
+        response.return_x = kwargs.pop('return_x', self.squeeze)
+
+        # Check for new labels
+        input_labels = kwargs.pop('input_labels', None)
+        if input_labels is not None:
+            response.input_labels = _process_labels(
+                input_labels, "input", response.ninputs)
+
+        output_labels = kwargs.pop('output_labels', None)
+        if output_labels is not None:
+            response.output_labels = _process_labels(
+                output_labels, "output", response.noutputs)
+
+        state_labels = kwargs.pop('state_labels', None)
+        if state_labels is not None:
+            response.state_labels = _process_labels(
+                state_labels, "state", response.nstates)
+
+        # Make sure no unknown keywords were passed
+        if len(kwargs) != 0:
+            raise ValueError("Unknown parameter(s) %s" % kwargs)
+
+        return response
+
+    @property
+    def time(self):
+
+        """Time vector.
+
+        Time values of the input/output response(s).
+
+        :type: 1D array"""
+        return self.t
+
+    # Getter for output (implements squeeze processing)
+    @property
+    def outputs(self):
+        """Time response output vector.
+
+        Output response of the system, indexed by either the output and time
+        (if only a single input is given) or the output, trace, and time
+        (for multiple traces).  See :attr:`TimeResponseData.squeeze` for a
+        description of how this can be modified using the `squeeze` keyword.
+
+        :type: 1D, 2D, or 3D array
+
+        """
+        t, y = _process_time_response(
+            self.t, self.y, issiso=self.issiso,
+            transpose=self.transpose, squeeze=self.squeeze)
+        return y
+
+    # Getter for states (implements squeeze processing)
+    @property
+    def states(self):
+        """Time response state vector.
+
+        Time evolution of the state vector, indexed indexed by either the
+        state and time (if only a single trace is given) or the state, trace,
+        and time (for multiple traces).  See :attr:`TimeResponseData.squeeze`
+        for a description of how this can be modified using the `squeeze`
+        keyword.
+
+        :type: 2D or 3D array
+
+        """
+        if self.x is None:
+            return None
+
+        elif self.squeeze is True:
+            x = self.x.squeeze()
+
+        elif self.ninputs == 1 and self.noutputs == 1 and \
+             self.ntraces == 1 and self.x.ndim == 3 and \
+             self.squeeze is not False:
+            # Single-input, single-output system with single trace
+            x = self.x[:, 0, :]
+
+        else:
+            # Return the full set of data
+            x = self.x
+
+        # Transpose processing
+        if self.transpose:
+            x = np.transpose(x, np.roll(range(x.ndim), 1))
+
+        return x
+
+    # Getter for inputs (implements squeeze processing)
+    @property
+    def inputs(self):
+        """Time response input vector.
+
+        Input(s) to the system, indexed by input (optiona), trace (optional),
+        and time.  If a 1D vector is passed, the input corresponds to a
+        scalar-valued input.  If a 2D vector is passed, then it can either
+        represent multiple single-input traces or a single multi-input trace.
+        The optional ``multi_trace`` keyword should be used to disambiguate
+        the two.  If a 3D vector is passed, then it represents a multi-trace,
+        multi-input signal, indexed by input, trace, and time.
+
+        See :attr:`TimeResponseData.squeeze` for a description of how the
+        dimensions of the input vector can be modified using the `squeeze`
+        keyword.
+
+        :type: 1D or 2D array
+
+        """
+        if self.u is None:
+            return None
+
+        t, u = _process_time_response(
+            self.t, self.u, issiso=self.issiso,
+            transpose=self.transpose, squeeze=self.squeeze)
+        return u
+
+    # Getter for legacy state (implements non-standard squeeze processing)
+    @property
+    def _legacy_states(self):
+        """Time response state vector (legacy version).
+
+        Time evolution of the state vector, indexed indexed by either the
+        state and time (if only a single trace is given) or the state,
+        trace, and time (for multiple traces).
+
+        The `legacy_states` property is not affected by the `squeeze` keyword
+        and hence it will always have these dimensions.
+
+        :type: 2D or 3D array
+
+        """
+
+        if self.x is None:
+            return None
+
+        elif self.ninputs == 1 and self.noutputs == 1 and \
+             self.ntraces == 1 and self.x.ndim == 3:
+            # Single-input, single-output system with single trace
+            x = self.x[:, 0, :]
+
+        else:
+            # Return the full set of data
+            x = self.x
+
+        # Transpose processing
+        if self.transpose:
+            x = np.transpose(x, np.roll(range(x.ndim), 1))
+
+        return x
+
+    # Implement iter to allow assigning to a tuple
+    def __iter__(self):
+        if not self.return_x:
+            return iter((self.time, self.outputs))
+        return iter((self.time, self.outputs, self._legacy_states))
+
+    # Implement (thin) getitem to allow access via legacy indexing
+    def __getitem__(self, index):
+        # See if we were passed a slice
+        if isinstance(index, slice):
+            if (index.start is None or index.start == 0) and index.stop == 2:
+                return (self.time, self.outputs)
+
+        # Otherwise assume we were passed a single index
+        if index == 0:
+            return self.time
+        if index == 1:
+            return self.outputs
+        if index == 2:
+            return self._legacy_states
+        raise IndexError
+
+    # Implement (thin) len to emulate legacy testing interface
+    def __len__(self):
+        return 3 if self.return_x else 2
+
+
+# Process signal labels
+def _process_labels(labels, signal, length):
+    """Process time response signal labels.
+
+    Parameters
+    ----------
+    labels : list of str or dict
+        Description of the labels for the signal.  This can be a list of
+        strings or a dict giving the index of each signal (used in iosys).
+
+    signal : str
+        Name of the signal being processed (for error messages).
+
+    length : int
+        Number of labels required.
+
+    Returns
+    -------
+    labels : list of str
+        List of labels.
+
+    """
+    if labels is None or len(labels) == 0:
+        return None
+
+    # See if we got passed a dictionary (from iosys)
+    if isinstance(labels, dict):
+        # Form inverse dictionary
+        ivd = {v: k for k, v in labels.items()}
+
+        try:
+            # Turn into a list
+            labels = [ivd[n] for n in range(len(labels))]
+        except KeyError:
+            raise ValueError("Name dictionary for %s is incomplete" % signal)
+
+    # Convert labels to a list
+    labels = list(labels)
+
+    # Make sure the signal list is the right length and type
+    if len(labels) != length:
+        raise ValueError("List of %s labels is the wrong length" % signal)
+    elif not all([isinstance(label, str) for label in labels]):
+        raise ValueError("List of %s labels must all  be strings" % signal)
+
+    return labels
 
 
 # Helper function for checking array-like parameters
 def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False,
                          transpose=False):
+
     """Helper function for checking array_like parameters.
 
     * Check type and shape of ``in_obj``.
@@ -237,9 +842,13 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         time simulations.
 
     return_x : bool, default=None
-        - If False, return only the time and output vectors.
-        - If True, also return the the state vector.
-        - If None, determine the returned variables by
+        Used if the time response data is assigned to a tuple:
+
+        * If False, return only the time and output vectors.
+
+        * If True, also return the the state vector.
+
+        * If None, determine the returned variables by
           config.defaults['forced_response.return_x'], which was True
           before version 0.9 and is False since then.
 
@@ -254,19 +863,25 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     Returns
     -------
-    T : array
-        Time values of the output.
+    results : TimeResponseData
+        Time response represented as a :class:`TimeResponseData` object
+        containing the following properties:
 
-    yout : array
-        Response of the system.  If the system is SISO and `squeeze` is not
-        True, the array is 1D (indexed by time).  If the system is not SISO or
-        `squeeze` is False, the array is 2D (indexed by the output number and
-        time).
+        * time (array): Time values of the output.
+
+        * outputs (array): Response of the system.  If the system is SISO and
+          `squeeze` is not True, the array is 1D (indexed by time).  If the
+          system is not SISO or `squeeze` is False, the array is 2D (indexed
+          by output and time).
 
-    xout : array
-        Time evolution of the state vector. Not affected by `squeeze`. Only
-        returned if `return_x` is True, or `return_x` is None and
-        config.defaults['forced_response.return_x'] is True.
+        * states (array): Time evolution of the state vector, represented as
+          a 2D array indexed by state and time.
+
+        * inputs (array): Input(s) to the system, indexed by input and time.
+
+        The return value of the system can also be accessed by assigning the
+        function to a tuple of length 2 (time, output) or of length 3 (time,
+        output, state) if ``return_x`` is ``True``.
 
     See Also
     --------
@@ -342,7 +957,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             # Make sure the input vector and time vector have same length
             if (U.ndim == 1 and U.shape[0] != T.shape[0]) or \
                     (U.ndim > 1 and U.shape[1] != T.shape[0]):
-                raise ValueError('Pamameter ``T`` must have same elements as'
+                raise ValueError('Parameter ``T`` must have same elements as'
                                  ' the number of columns in input array ``U``')
             if U.ndim == 0:
                 U = np.full((n_inputs, T.shape[0]), U)
@@ -368,6 +983,13 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
                               'Parameter ``X0``: ', squeeze=True)
 
+    # Test if U has correct shape and type
+    legal_shapes = [(n_steps,), (1, n_steps)] if n_inputs == 1 else \
+        [(n_inputs, n_steps)]
+    U = _check_convert_array(U, legal_shapes,
+                             'Parameter ``U``: ', squeeze=False,
+                             transpose=transpose)
+
     xout = np.zeros((n_states, n_steps))
     xout[:, 0] = X0
     yout = np.zeros((n_outputs, n_steps))
@@ -388,17 +1010,11 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
         # General algorithm that interpolates U in between output points
         else:
-            # Test if U has correct shape and type
-            legal_shapes = [(n_steps,), (1, n_steps)] if n_inputs == 1 else \
-                           [(n_inputs, n_steps)]
-            U = _check_convert_array(U, legal_shapes,
-                                     'Parameter ``U``: ', squeeze=False,
-                                     transpose=transpose)
-            # convert 1D array to 2D array with only one row
-            if len(U.shape) == 1:
+            # convert input from 1D array to 2D array with only one row
+            if U.ndim == 1:
                 U = U.reshape(1, -1)  # pylint: disable=E1103
 
-            # Algorithm: to integrate from time 0 to time dt, with linear
+        # Algorithm: to integrate from time 0 to time dt, with linear
             # interpolation between inputs u(0) = u0 and u(dt) = u1, we solve
             #   xdot = A x + B u,        x(0) = x0
             #   udot = (u1 - u0) / dt,   u(0) = u0.
@@ -471,29 +1087,29 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             tout = T            # Return exact list of time steps
             yout = yout[::inc, :]
             xout = xout[::inc, :]
+        else:
+            # Interpolate the input to get the right number of points
+            U = sp.interpolate.interp1d(T, U)(tout)
 
         # Transpose the output and state vectors to match local convention
         xout = np.transpose(xout)
         yout = np.transpose(yout)
 
-    return _process_time_response(sys, tout, yout, xout, transpose=transpose,
-                                  return_x=return_x, squeeze=squeeze)
+    return TimeResponseData(
+        tout, yout, xout, U, issiso=sys.issiso(),
+        transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
 # Process time responses in a uniform way
 def _process_time_response(
-        sys, tout, yout, xout, transpose=None, return_x=False,
-        squeeze=None, input=None, output=None):
+        tout, yout, issiso=False, transpose=None, squeeze=None):
     """Process time response signals.
 
-    This function processes the outputs of the time response functions and
-    processes the transpose and squeeze keywords.
+    This function processes the outputs (or inputs) of time response
+    functions and processes the transpose and squeeze keywords.
 
     Parameters
     ----------
-    sys : LTI or InputOutputSystem
-        System that generated the data (used to check if SISO/MIMO).
-
     T : 1D array
         Time values of the output.  Ignored if None.
 
@@ -503,20 +1119,15 @@ def _process_time_response(
         systems with no input indexing, such as initial_response or forced
         response) or a 3D array indexed by output, input, and time.
 
-    xout : array, optional
-        Individual response of each x variable (if return_x is True). For a
-        SISO system (or if a single input is specified), this should be a 2D
-        array indexed by the state index and time (for single input systems)
-        or a 3D array indexed by state, input, and time. Ignored if None.
+    issiso : bool, optional
+        If ``True``, process data as single-input, single-output data.
+        Default is ``False``.
 
     transpose : bool, optional
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`).  Default
         value is False.
 
-    return_x : bool, optional
-        If True, return the state vector (default = False).
-
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then the
         output response is returned as a 1D array (indexed by time).  If
@@ -526,16 +1137,10 @@ def _process_time_response(
         the system is SISO. The default value can be set using
         config.defaults['control.squeeze_time_response'].
 
-    input : int, optional
-        If present, the response represents only the listed input.
-
-    output : int, optional
-        If present, the response represents only the listed output.
-
     Returns
     -------
     T : 1D array
-        Time values of the output
+        Time values of the output.
 
     yout : ndarray
         Response of the system.  If the system is SISO and squeeze is not
@@ -543,20 +1148,11 @@ def _process_time_response(
         squeeze is False, the array is either 2D (indexed by output and time)
         or 3D (indexed by input, output, and time).
 
-    xout : array, optional
-        Individual response of each x variable (if return_x is True). For a
-        SISO system (or if a single input is specified), xout is a 2D array
-        indexed by the state index and time.  For a non-SISO system, xout is a
-        3D array indexed by the state, the input, and time.  The shape of xout
-        is not affected by the ``squeeze`` keyword.
     """
     # If squeeze was not specified, figure out the default (might remain None)
     if squeeze is None:
         squeeze = config.defaults['control.squeeze_time_response']
 
-    # Determine if the system is SISO
-    issiso = sys.issiso() or (input is not None and output is not None)
-
     # Figure out whether and how to squeeze output data
     if squeeze is True:         # squeeze all dimensions
         yout = np.squeeze(yout)
@@ -564,16 +1160,12 @@ def _process_time_response(
         pass
     elif squeeze is None:       # squeeze signals if SISO
         if issiso:
-            if len(yout.shape) == 3:
+            if yout.ndim == 3:
                 yout = yout[0][0]       # remove input and output
             else:
                 yout = yout[0]          # remove input
     else:
-        raise ValueError("unknown squeeze value")
-
-    # Figure out whether and how to squeeze the state data
-    if issiso and xout is not None and len(xout.shape) > 2:
-        xout = xout[:, 0, :]            # remove input
+        raise ValueError("Unknown squeeze value")
 
     # See if we need to transpose the data back into MATLAB form
     if transpose:
@@ -582,11 +1174,9 @@ def _process_time_response(
 
         # For signals, put the last index (time) into the first slot
         yout = np.transpose(yout, np.roll(range(yout.ndim), 1))
-        if xout is not None:
-            xout = np.transpose(xout, np.roll(range(xout.ndim), 1))
 
     # Return time, output, and (optionally) state
-    return (tout, yout, xout) if return_x else (tout, yout)
+    return tout, yout
 
 
 def _get_ss_simo(sys, input=None, output=None, squeeze=None):
@@ -609,7 +1199,7 @@ def _get_ss_simo(sys, input=None, output=None, squeeze=None):
     sys_ss = _convert_to_statespace(sys)
     if sys_ss.issiso():
         return squeeze, sys_ss
-    elif squeeze == None and (input is None or output is None):
+    elif squeeze is None and (input is None or output is None):
         # Don't squeeze outputs if resulting system turns out to be siso
         # Note: if we expand input to allow a tuple, need to update this check
         squeeze = False
@@ -662,7 +1252,8 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
 
     input : int, optional
         Only compute the step response for the listed input.  If not
-        specified, the step responses for each independent input are computed.
+        specified, the step responses for each independent input are
+        computed (as separate traces).
 
     output : int, optional
         Only report the step response for the listed output.  If not
@@ -678,7 +1269,8 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         value is False.
 
     return_x : bool, optional
-        If True, return the state vector (default = False).
+        If True, return the state vector when assigning to a tuple (default =
+        False).  See :func:`forced_response` for more details.
 
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then the
@@ -691,21 +1283,27 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
 
     Returns
     -------
-    T : 1D array
-        Time values of the output
+    results : TimeResponseData
+        Time response represented as a :class:`TimeResponseData` object
+        containing the following properties:
 
-    yout : ndarray
-        Response of the system.  If the system is SISO and squeeze is not
-        True, the array is 1D (indexed by time).  If the system is not SISO or
-        squeeze is False, the array is 3D (indexed by the input, output, and
-        time).
+        * time (array): Time values of the output.
+
+        * outputs (array): Response of the system.  If the system is SISO and
+          squeeze is not True, the array is 1D (indexed by time).  If the
+          system is not SISO or ``squeeze`` is False, the array is 3D (indexed
+          by the output, trace, and time).
+
+        * states (array): Time evolution of the state vector, represented as
+          either a 2D array indexed by state and time (if SISO) or a 3D array
+          indexed by state, trace, and time.  Not affected by ``squeeze``.
 
-    xout : array, optional
-        Individual response of each x variable (if return_x is True). For a
-        SISO system (or if a single input is specified), xout is a 2D array
-        indexed by the state index and time.  For a non-SISO system, xout is a
-        3D array indexed by the state, the input, and time.  The shape of xout
-        is not affected by the ``squeeze`` keyword.
+        * inputs (array): Input(s) to the system, indexed in the same manner
+          as ``outputs``.
+
+        The return value of the system can also be accessed by assigning the
+        function to a tuple of length 2 (time, output) or of length 3 (time,
+        output, state) if ``return_x`` is ``True``.
 
     See Also
     --------
@@ -741,6 +1339,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
     noutputs = sys.noutputs if output is None else 1
     yout = np.empty((noutputs, ninputs, np.asarray(T).size))
     xout = np.empty((sys.nstates, ninputs, np.asarray(T).size))
+    uout = np.empty((ninputs, ninputs, np.asarray(T).size))
 
     # Simulate the response for each input
     for i in range(sys.ninputs):
@@ -751,16 +1350,18 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         # Create a set of single inputs system for simulation
         squeeze, simo = _get_ss_simo(sys, i, output, squeeze=squeeze)
 
-        out = forced_response(simo, T, U, X0, transpose=False,
-                              return_x=return_x, squeeze=True)
+        response = forced_response(simo, T, U, X0, squeeze=True)
         inpidx = i if input is None else 0
-        yout[:, inpidx, :] = out[1]
-        if return_x:
-            xout[:, i, :] = out[2]
+        yout[:, inpidx, :] = response.y
+        xout[:, inpidx, :] = response.x
+        uout[:, inpidx, :] = U
+
+    # Figure out if the system is SISO or not
+    issiso = sys.issiso() or (input is not None and output is not None)
 
-    return _process_time_response(
-        sys, out[0], yout, xout, transpose=transpose, return_x=return_x,
-        squeeze=squeeze, input=input, output=output)
+    return TimeResponseData(
+        response.time, yout, xout, uout, issiso=issiso,
+        transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
 def step_info(sysdata, T=None, T_num=None, yfinal=None,
@@ -985,6 +1586,7 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None,
 
     return ret[0][0] if retsiso else ret
 
+
 def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
                      transpose=False, return_x=False, squeeze=None):
     # pylint: disable=W0622
@@ -1028,7 +1630,8 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         value is False.
 
     return_x : bool, optional
-        If True, return the state vector (default = False).
+        If True, return the state vector when assigning to a tuple (default =
+        False).  See :func:`forced_response` for more details.
 
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then the
@@ -1041,17 +1644,24 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
 
     Returns
     -------
-    T : array
-        Time values of the output
+    results : TimeResponseData
+        Time response represented as a :class:`TimeResponseData` object
+        containing the following properties:
 
-    yout : array
-        Response of the system.  If the system is SISO and squeeze is not
-        True, the array is 1D (indexed by time).  If the system is not SISO or
-        squeeze is False, the array is 2D (indexed by the output number and
-        time).
+        * time (array): Time values of the output.
+
+        * outputs (array): Response of the system.  If the system is SISO and
+          squeeze is not True, the array is 1D (indexed by time).  If the
+          system is not SISO or ``squeeze`` is False, the array is 2D (indexed
+          by the output and time).
+
+        * states (array): Time evolution of the state vector, represented as
+          either a 2D array indexed by state and time (if SISO).  Not affected
+          by ``squeeze``.
 
-    xout : array, optional
-        Individual response of each x variable (if return_x is True).
+        The return value of the system can also be accessed by assigning the
+        function to a tuple of length 2 (time, output) or of length 3 (time,
+        output, state) if ``return_x`` is ``True``.
 
     See Also
     --------
@@ -1073,10 +1683,17 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
     # The initial vector X0 is created in forced_response(...) if necessary
     if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)
-    U = np.zeros_like(T)
 
-    return forced_response(sys, T, U, X0, transpose=transpose,
-                           return_x=return_x, squeeze=squeeze)
+    # Compute the forced response
+    response = forced_response(sys, T, 0, X0)
+
+    # Figure out if the system is SISO or not
+    issiso = sys.issiso() or (input is not None and output is not None)
+
+    # Store the response without an input
+    return TimeResponseData(
+        response.t, response.y, response.x, None, issiso=issiso,
+        transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
 def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
@@ -1126,7 +1743,8 @@ def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         value is False.
 
     return_x : bool, optional
-        If True, return the state vector (default = False).
+        If True, return the state vector when assigning to a tuple (default =
+        False).  See :func:`forced_response` for more details.
 
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then the
@@ -1139,21 +1757,24 @@ def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
 
     Returns
     -------
-    T : array
-        Time values of the output
+    results : TimeResponseData
+        Impulse response represented as a :class:`TimeResponseData` object
+        containing the following properties:
 
-    yout : array
-        Response of the system.  If the system is SISO and squeeze is not
-        True, the array is 1D (indexed by time).  If the system is not SISO or
-        squeeze is False, the array is 2D (indexed by the output number and
-        time).
+        * time (array): Time values of the output.
+
+        * outputs (array): Response of the system.  If the system is SISO and
+          squeeze is not True, the array is 1D (indexed by time).  If the
+          system is not SISO or ``squeeze`` is False, the array is 3D (indexed
+          by the output, trace, and time).
 
-    xout : array, optional
-        Individual response of each x variable (if return_x is True). For a
-        SISO system (or if a single input is specified), xout is a 2D array
-        indexed by the state index and time.  For a non-SISO system, xout is a
-        3D array indexed by the state, the input, and time.  The shape of xout
-        is not affected by the ``squeeze`` keyword.
+        * states (array): Time evolution of the state vector, represented as
+          either a 2D array indexed by state and time (if SISO) or a 3D array
+          indexed by state, trace, and time.  Not affected by ``squeeze``.
+
+        The return value of the system can also be accessed by assigning the
+        function to a tuple of length 2 (time, output) or of length 3 (time,
+        output, state) if ``return_x`` is ``True``.
 
     See Also
     --------
@@ -1196,6 +1817,7 @@ def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
     noutputs = sys.noutputs if output is None else 1
     yout = np.empty((noutputs, ninputs, np.asarray(T).size))
     xout = np.empty((sys.nstates, ninputs, np.asarray(T).size))
+    uout = np.full((ninputs, ninputs, np.asarray(T).size), None)
 
     # Simulate the response for each input
     for i in range(sys.ninputs):
@@ -1218,21 +1840,22 @@ def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
             new_X0 = B + X0
         else:
             new_X0 = X0
-            U[0] = 1./simo.dt # unit area impulse
+            U[0] = 1./simo.dt           # unit area impulse
 
         # Simulate the impulse response fo this input
-        out = forced_response(simo, T, U, new_X0, transpose=False,
-                              return_x=return_x, squeeze=squeeze)
+        response = forced_response(simo, T, U, new_X0)
 
         # Store the output (and states)
         inpidx = i if input is None else 0
-        yout[:, inpidx, :] = out[1]
-        if return_x:
-            xout[:, i, :] = out[2]
+        yout[:, inpidx, :] = response.y
+        xout[:, inpidx, :] = response.x
+
+    # Figure out if the system is SISO or not
+    issiso = sys.issiso() or (input is not None and output is not None)
 
-    return _process_time_response(
-        sys, out[0], yout, xout, transpose=transpose, return_x=return_x,
-        squeeze=squeeze, input=input, output=output)
+    return TimeResponseData(
+        response.time, yout, xout, uout, issiso=issiso,
+        transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
 # utility function to find time period and time increment using pole locations
@@ -1283,11 +1906,11 @@ def _ideal_tfinal_and_dt(sys, is_step=True):
     """
 
     sqrt_eps = np.sqrt(np.spacing(1.))
-    default_tfinal = 5  # Default simulation horizon
+    default_tfinal = 5                  # Default simulation horizon
     default_dt = 0.1
-    total_cycles = 5  # number of cycles for oscillating modes
-    pts_per_cycle = 25  # Number of points divide a period of oscillation
-    log_decay_percent = np.log(1000)  # Factor of reduction for real pole decays
+    total_cycles = 5                    # Number cycles for oscillating modes
+    pts_per_cycle = 25                  # Number points divide period of osc
+    log_decay_percent = np.log(1000)    # Reduction factor for real pole decays
 
     if sys._isstatic():
         tfinal = default_tfinal
@@ -1333,13 +1956,15 @@ def _ideal_tfinal_and_dt(sys, is_step=True):
 
         if p_int.size > 0:
             tfinal = tfinal * 5
-    else: # cont time
+    else:       # cont time
         sys_ss = _convert_to_statespace(sys)
         # Improve conditioning via balancing and zeroing tiny entries
-        # See <w,v> for [[1,2,0], [9,1,0.01], [1,2,10*np.pi]] before/after balance
+        # See <w,v> for [[1,2,0], [9,1,0.01], [1,2,10*np.pi]]
+        #   before/after balance
         b, (sca, perm) = matrix_balance(sys_ss.A, separate=True)
         p, l, r = eig(b, left=True, right=True)
-        # Reciprocal of inner product <w,v> for each eigval, (bound the ~infs by 1e12)
+        # Reciprocal of inner product <w,v> for each eigval, (bound the
+        #   ~infs by 1e12)
         # G = Transfer([1], [1,0,1]) gives zero sensitivity (bound by 1e-12)
         eig_sens = np.reciprocal(maximum(1e-12, einsum('ij,ij->j', l, r).real))
         eig_sens = minimum(1e12, eig_sens)
@@ -1359,9 +1984,9 @@ def _ideal_tfinal_and_dt(sys, is_step=True):
         dc = np.zeros_like(p, dtype=float)
         # well-conditioned nonzero poles, np.abs just in case
         ok = np.abs(eig_sens) <= 1/sqrt_eps
-        # the averaged t->inf response of each simple eigval on each i/o channel
-        # See, A = [[-1, k], [0, -2]], response sizes are k-dependent (that is
-        # R/L eigenvector dependent)
+        # the averaged t->inf response of each simple eigval on each i/o
+        # channel. See, A = [[-1, k], [0, -2]], response sizes are
+        # k-dependent (that is R/L eigenvector dependent)
         dc[ok] = norm(v[ok, :], axis=1)*norm(w[:, ok], axis=0)*eig_sens[ok]
         dc[wn != 0.] /= wn[wn != 0] if is_step else 1.
         dc[wn == 0.] = 0.
@@ -1384,8 +2009,10 @@ def _ideal_tfinal_and_dt(sys, is_step=True):
         # The rest ~ts = log(%ss value) / exp(Re(eigval)t)
         texp_mode = log_decay_percent / np.abs(psub[~iw & ~ints].real)
         tfinal += texp_mode.tolist()
-        dt += minimum(texp_mode / 50,
-                    (2 * np.pi / pts_per_cycle / wnsub[~iw & ~ints])).tolist()
+        dt += minimum(
+            texp_mode / 50,
+            (2 * np.pi / pts_per_cycle / wnsub[~iw & ~ints])
+        ).tolist()
 
         # All integrators?
         if len(tfinal) == 0:
@@ -1396,13 +2023,14 @@ def _ideal_tfinal_and_dt(sys, is_step=True):
 
     return tfinal, dt
 
+
 def _default_time_vector(sys, N=None, tfinal=None, is_step=True):
     """Returns a time vector that has a reasonable number of points.
     if system is discrete-time, N is ignored """
 
     N_max = 5000
-    N_min_ct = 100 # min points for cont time systems
-    N_min_dt = 20 # more common to see just a few samples in discrete-time
+    N_min_ct = 100    # min points for cont time systems
+    N_min_dt = 20     # more common to see just a few samples in discrete time
 
     ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys, is_step=is_step)
 
@@ -1415,7 +2043,7 @@ def _default_time_vector(sys, N=None, tfinal=None, is_step=True):
             tfinal = sys.dt * (N-1)
         else:
             N = int(np.ceil(tfinal/sys.dt)) + 1
-            tfinal = sys.dt * (N-1) # make tfinal an integer multiple of sys.dt
+            tfinal = sys.dt * (N-1)  # make tfinal integer multiple of sys.dt
     else:
         if tfinal is None:
             # for continuous time, simulate to ideal_tfinal but limit N
diff --git a/doc/classes.rst b/doc/classes.rst
index b80b7dd54..0753271c4 100644
--- a/doc/classes.rst
+++ b/doc/classes.rst
@@ -17,6 +17,7 @@ these directly.
    TransferFunction
    StateSpace
    FrequencyResponseData
+   TimeResponseData
 
 Input/output system subclasses
 ==============================
diff --git a/doc/conventions.rst b/doc/conventions.rst
index 63f3fac2c..462a71408 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -134,13 +134,12 @@ Types:
     * **Arguments** can be **arrays**, **matrices**, or **nested lists**.
     * **Return values** are **arrays** (not matrices).
 
-The time vector is either 1D, or 2D with shape (1, n)::
+The time vector is a 1D array with shape (n, )::
 
-      T = [[t1,     t2,     t3,     ..., tn    ]]
+      T = [t1,     t2,     t3,     ..., tn    ]
 
 Input, state, and output all follow the same convention. Columns are different
-points in time, rows are different components. When there is only one row, a
-1D object is accepted or returned, which adds convenience for SISO systems::
+points in time, rows are different components::
 
       U = [[u1(t1), u1(t2), u1(t3), ..., u1(tn)]
            [u2(t1), u2(t2), u2(t3), ..., u2(tn)]
@@ -153,6 +152,9 @@ points in time, rows are different components. When there is only one row, a
 So, U[:,2] is the system's input at the third point in time; and U[1] or U[1,:]
 is the sequence of values for the system's second input.
 
+When there is only one row, a 1D object is accepted or returned, which adds
+convenience for SISO systems:
+
 The initial conditions are either 1D, or 2D with shape (j, 1)::
 
      X0 = [[x1]
@@ -161,23 +163,43 @@ The initial conditions are either 1D, or 2D with shape (j, 1)::
            ...
            [xj]]
 
-As all simulation functions return *arrays*, plotting is convenient::
+Functions that return time responses (e.g., :func:`forced_response`,
+:func:`impulse_response`, :func:`input_output_response`,
+:func:`initial_response`, and :func:`step_response`) return a
+:class:`TimeResponseData` object that contains the data for the time
+response.  These data can be accessed via the ``time``, ``outputs``,
+``states`` and ``inputs`` properties::
+
+    sys = rss(4, 1, 1)
+    response = step_response(sys)
+    plot(response.time, response.outputs)
+
+The dimensions of the response properties depend on the function being
+called and whether the system is SISO or MIMO.  In addition, some time
+response function can return multiple "traces" (input/output pairs),
+such as the :func:`step_response` function applied to a MIMO system,
+which will compute the step response for each input/output pair.  See
+:class:`TimeResponseData` for more details.
+
+The time response functions can also be assigned to a tuple, which extracts
+the time and output (and optionally the state, if the `return_x` keyword is
+used).  This allows simple commands for plotting::
 
     t, y = step_response(sys)
     plot(t, y)
 
 The output of a MIMO system can be plotted like this::
 
-    t, y = forced_response(sys, u, t)
+    t, y = forced_response(sys, t, u)
     plot(t, y[0], label='y_0')
     plot(t, y[1], label='y_1')
 
-The convention also works well with the state space form of linear systems. If
-``D`` is the feedthrough *matrix* of a linear system, and ``U`` is its input
-(*matrix* or *array*), then the feedthrough part of the system's response,
-can be computed like this::
+The convention also works well with the state space form of linear
+systems. If ``D`` is the feedthrough matrix (2D array) of a linear system,
+and ``U`` is its input (array), then the feedthrough part of the system's
+response, can be computed like this::
 
-    ft = D * U
+    ft = D @ U
 
 
 .. currentmodule:: control
@@ -210,27 +232,29 @@ on standard configurations.
 
 Selected variables that can be configured, along with their default values:
 
-  * freqplot.dB (False): Bode plot magnitude plotted in dB (otherwise powers of 10)
+  * freqplot.dB (False): Bode plot magnitude plotted in dB (otherwise powers
+    of 10)
     
   * freqplot.deg (True): Bode plot phase plotted in degrees (otherwise radians)
     
-  * freqplot.Hz (False): Bode plot frequency plotted in Hertz (otherwise rad/sec)
+  * freqplot.Hz (False): Bode plot frequency plotted in Hertz (otherwise
+    rad/sec)
     
   * freqplot.grid (True): Include grids for magnitude and phase plots
     
   * freqplot.number_of_samples (1000): Number of frequency points in Bode plots
     
-  * freqplot.feature_periphery_decade (1.0): How many decades to include in the
-    frequency range on both sides of features (poles, zeros).
+  * freqplot.feature_periphery_decade (1.0): How many decades to include in
+    the frequency range on both sides of features (poles, zeros).
 
-  * statesp.use_numpy_matrix (True): set the return type for state space matrices to
-    `numpy.matrix` (verus numpy.ndarray)
+  * statesp.use_numpy_matrix (True): set the return type for state space
+    matrices to `numpy.matrix` (verus numpy.ndarray)
 
-  * statesp.default_dt and xferfcn.default_dt (None): set the default value of dt when
-    constructing new LTI systems
+  * statesp.default_dt and xferfcn.default_dt (None): set the default value
+    of dt when constructing new LTI systems
 
-  * statesp.remove_useless_states (True): remove states that have no effect on the
-    input-output dynamics of the system
+  * statesp.remove_useless_states (True): remove states that have no effect
+    on the input-output dynamics of the system
 
 Additional parameter variables are documented in individual functions