diff --git a/control/__init__.py b/control/__init__.py
index 9781ab80e..120d16325 100644
--- a/control/__init__.py
+++ b/control/__init__.py
@@ -77,6 +77,10 @@
 from .xferfcn import *
 from .frdata import *
 
+# Time responses and plotting
+from .timeresp import *
+from .timeplot import *
+
 from .bdalg import *
 from .delay import *
 from .descfcn import *
@@ -91,7 +95,6 @@
 from .rlocus import *
 from .statefbk import *
 from .stochsys import *
-from .timeresp import *
 from .ctrlutil import *
 from .canonical import *
 from .robust import *
diff --git a/control/config.py b/control/config.py
index 987693c2d..1ed8b5dd5 100644
--- a/control/config.py
+++ b/control/config.py
@@ -135,6 +135,9 @@ def reset_defaults():
     from .optimal import _optimal_defaults
     defaults.update(_optimal_defaults)
 
+    from .timeplot import _timeplot_defaults
+    defaults.update(_timeplot_defaults)
+
 
 def _get_param(module, param, argval=None, defval=None, pop=False, last=False):
     """Return the default value for a configuration option.
diff --git a/control/nlsys.py b/control/nlsys.py
index c540decb6..f11d08a31 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -1488,6 +1488,7 @@ def ufun(t):
         return TimeResponseData(
             t_eval, y, None, u, issiso=sys.issiso(),
             output_labels=sys.output_labels, input_labels=sys.input_labels,
+            title="Input/output response for " + sys.name, sysname=sys.name,
             transpose=transpose, return_x=return_x, squeeze=squeeze)
 
     # Create a lambda function for the right hand side
@@ -1567,7 +1568,8 @@ def ivp_rhs(t, x):
     return TimeResponseData(
         soln.t, y, soln.y, u, issiso=sys.issiso(),
         output_labels=sys.output_labels, input_labels=sys.input_labels,
-        state_labels=sys.state_labels,
+        state_labels=sys.state_labels, sysname=sys.name,
+        title="Input/output response for " + sys.name,
         transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index 1abc7547e..19f4bb627 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -26,6 +26,7 @@
 import control.tests.statefbk_test as statefbk_test
 import control.tests.stochsys_test as stochsys_test
 import control.tests.trdata_test as trdata_test
+import control.tests.timeplot_test as timeplot_test
 
 @pytest.mark.parametrize("module, prefix", [
     (control, ""), (control.flatsys, "flatsys."), (control.optimal, "optimal.")
@@ -74,7 +75,8 @@ def test_kwarg_search(module, prefix):
                 # @parametrize messes up the check, but we know it is there
                 pass
 
-            elif source and source.find('unrecognized keyword') < 0:
+            elif source and source.find('unrecognized keyword') < 0 and \
+                 source.find('unexpected keyword') < 0:
                 warnings.warn(
                     f"'unrecognized keyword' not found in unit test "
                     f"for {name}")
@@ -161,7 +163,21 @@ def test_matplotlib_kwargs(function, nsysargs, moreargs, kwargs, mplcleanup):
         function(*args, **kwargs, unknown=None)
 
 
+@pytest.mark.parametrize(
+    "function", [control.time_response_plot, control.TimeResponseData.plot])
+def test_time_response_plot_kwargs(function):
+    # Create a system for testing
+    response = control.step_response(control.rss(4, 2, 2))
+
+    # Call the plotting function normally and make sure it works
+    function(response)
 
+    # Now add an unrecognized keyword and make sure there is an error
+    with pytest.raises(AttributeError,
+                       match="(has no property|unexpected keyword)"):
+        function(response, unknown=None)
+    
+    
 #
 # List of all unit tests that check for unrecognized keywords
 #
@@ -185,6 +201,7 @@ def test_matplotlib_kwargs(function, nsysargs, moreargs, kwargs, mplcleanup):
     'gangof4_plot': test_matplotlib_kwargs,
     'input_output_response': test_unrecognized_kwargs,
     'interconnect': interconnect_test.test_interconnect_exceptions,
+    'time_response_plot': timeplot_test.test_errors,
     'linearize': test_unrecognized_kwargs,
     'lqe': test_unrecognized_kwargs,
     'lqr': test_unrecognized_kwargs,
@@ -230,6 +247,7 @@ def test_matplotlib_kwargs(function, nsysargs, moreargs, kwargs, mplcleanup):
     'StateSpace.__init__': test_unrecognized_kwargs,
     'StateSpace.sample': test_unrecognized_kwargs,
     'TimeResponseData.__call__': trdata_test.test_response_copy,
+    'TimeResponseData.plot': timeplot_test.test_errors,
     'TransferFunction.__init__': test_unrecognized_kwargs,
     'TransferFunction.sample': test_unrecognized_kwargs,
     'optimal.OptimalControlProblem.__init__':
diff --git a/control/tests/timeplot_test.py b/control/tests/timeplot_test.py
new file mode 100644
index 000000000..7cdde5c54
--- /dev/null
+++ b/control/tests/timeplot_test.py
@@ -0,0 +1,511 @@
+# timeplot_test.py - test out time response plots
+# RMM, 23 Jun 2023
+
+import pytest
+import control as ct
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+from control.tests.conftest import slycotonly
+
+# Detailed test of (almost) all functionality
+#
+# The commented out rows lead to very long testing times => these should be
+# used only for developmental testing and not day-to-day testing.
+@pytest.mark.parametrize(
+    "sys", [
+        # ct.rss(1, 1, 1, strictly_proper=True, name="rss"),
+        ct.nlsys(
+            lambda t, x, u, params: -x + u, None,
+            inputs=1, outputs=1, states=1, name="nlsys"),
+        # ct.rss(2, 1, 2, strictly_proper=True, name="rss"),
+        ct.rss(2, 2, 1, strictly_proper=True, name="rss"),
+        # ct.drss(2, 2, 2, name="drss"),
+        # ct.rss(2, 2, 3, strictly_proper=True, name="rss"),
+    ])
+# @pytest.mark.parametrize("transpose", [False, True])
+# @pytest.mark.parametrize("plot_inputs", [False, None, True, 'overlay'])
+# @pytest.mark.parametrize("plot_outputs", [True, False])
+# @pytest.mark.parametrize("overlay_signals", [False, True])
+# @pytest.mark.parametrize("overlay_traces", [False, True])
+# @pytest.mark.parametrize("second_system", [False, True])
+# @pytest.mark.parametrize("fcn", [
+#     ct.step_response, ct.impulse_response, ct.initial_response,
+#     ct.forced_response])
+@pytest.mark.parametrize(       # combinatorial-style test (faster)
+    "fcn,                  pltinp,    pltout, cmbsig, cmbtrc, trpose, secsys",
+    [(ct.step_response,    False,     True,   False,  False,  False,  False),
+     (ct.step_response,    None,      True,   False,  False,  False,  False),
+     (ct.step_response,    True,      True,   False,  False,  False,  False),
+     (ct.step_response,    'overlay', True,   False,  False,  False,  False),
+     (ct.step_response,    'overlay', True,   True,   False,  False,  False),
+     (ct.step_response,    'overlay', True,   False,  True,   False,  False),
+     (ct.step_response,    'overlay', True,   False,  False,  True,   False),
+     (ct.step_response,    'overlay', True,   False,  False,  False,  True),
+     (ct.step_response,    False,     False,  False,  False,  False,  False),
+     (ct.step_response,    None,      False,  False,  False,  False,  False),
+     (ct.step_response,    'overlay', False,  False,  False,  False,  False),
+     (ct.step_response,    True,      True,   False,  True,   False,  False),
+     (ct.step_response,    True,      True,   False,  False,  False,  True),
+     (ct.step_response,    True,      True,   False,  True,   False,  True),
+     (ct.step_response,    True,      True,   True,   False,   True,   True),
+     (ct.step_response,    True,      True,   False,  True,    True,   True),
+     (ct.impulse_response, False,     True,   True,   False,  False,  False),
+     (ct.initial_response, None,      True,   False,  False,  False,  False),
+     (ct.initial_response, False,     True,   False,  False,  False,  False),
+     (ct.initial_response, True,      True,   False,  False,  False,  False),
+     (ct.forced_response,  True,      True,   False,  False,  False,  False),
+     (ct.forced_response,  None,      True,   False,  False,  False,  False),
+     (ct.forced_response,  False,     True,   False,  False,  False,  False),
+     (ct.forced_response,  True,      True,   True,   False,  False,  False),
+     (ct.forced_response,  True,      True,   True,   True,   False,  False),
+     (ct.forced_response,  True,      True,   True,   True,   True,   False),
+     (ct.forced_response,  True,      True,   True,   True,   True,   True),
+     (ct.forced_response,  'overlay', True,   True,   True,   False,  True),
+     (ct.input_output_response,
+                           True,      True,   False,  False,  False,  False),
+     ])
+
+def test_response_plots(
+        fcn, sys, pltinp, pltout, cmbsig, cmbtrc,
+        trpose, secsys, clear=True):
+    # Figure out the time range to use and check some special cases
+    if not isinstance(sys, ct.lti.LTI):
+        if fcn == ct.impulse_response:
+            pytest.skip("impulse response not implemented for nlsys")
+
+        # Nonlinear systems require explicit time limits
+        T = 10
+        timepts = np.linspace(0, T)
+
+    elif isinstance(sys, ct.TransferFunction) and fcn == ct.initial_response:
+        pytest.skip("initial response not tested for tf")
+
+    else:
+        # Linear systems figure things out on their own
+        T = None
+        timepts = np.linspace(0, 10)    # for input_output_response
+
+    # Save up the keyword arguments
+    kwargs = dict(
+        plot_inputs=pltinp, plot_outputs=pltout, transpose=trpose,
+        overlay_signals=cmbsig, overlay_traces=cmbtrc)
+
+    # Create the response
+    if fcn is ct.input_output_response and \
+       not isinstance(sys, ct.NonlinearIOSystem):
+        # Skip transfer functions and other non-state space systems
+        return None
+    if fcn in [ct.input_output_response, ct.forced_response]:
+        U = np.zeros((sys.ninputs, timepts.size))
+        for i in range(sys.ninputs):
+            U[i] = np.cos(timepts * i + i)
+        args = [timepts, U]
+
+    elif fcn == ct.initial_response:
+        args = [T, np.ones(sys.nstates)]   # T, X0
+
+    elif not isinstance(sys, ct.lti.LTI):
+        args = [T]              # nonlinear systems require final time
+
+    else:                       # step, initial, impulse responses
+        args = []
+
+    # Create a new figure (in case previous one is of the same size) and plot
+    if not clear:
+        plt.figure()
+    response = fcn(sys, *args)
+
+    # Look for cases where there are no data to plot
+    if not pltout and (
+            pltinp is False or response.ninputs == 0 or
+            pltinp is None and response.plot_inputs is False):
+        with pytest.raises(ValueError, match=".* no data to plot"):
+            out = response.plot(**kwargs)
+        return None
+    elif not pltout and pltinp == 'overlay':
+        with pytest.raises(ValueError, match="can't overlay inputs"):
+            out = response.plot(**kwargs)
+        return None
+    elif pltinp in [True, 'overlay'] and response.ninputs == 0:
+        with pytest.raises(ValueError, match=".* but no inputs"):
+            out = response.plot(**kwargs)
+        return None
+
+    out = response.plot(**kwargs)
+
+    # Make sure all of the outputs are of the right type
+    nlines_plotted = 0
+    for ax_lines in np.nditer(out, flags=["refs_ok"]):
+        for line in ax_lines.item():
+            assert isinstance(line, mpl.lines.Line2D)
+            nlines_plotted += 1
+
+    # Make sure number of plots is correct
+    if pltinp is None:
+        if fcn in [ct.forced_response, ct.input_output_response]:
+            pltinp = True
+        else:
+            pltinp = False
+    ntraces = max(1, response.ntraces)
+    nlines_expected = (response.ninputs if pltinp else 0) * ntraces + \
+        (response.noutputs if pltout else 0) * ntraces
+    assert nlines_plotted == nlines_expected
+
+    # Save the old axes to compare later
+    old_axes = plt.gcf().get_axes()
+
+    # Add additional data (and provide info in the title)
+    if secsys:
+        newsys = ct.rss(
+            sys.nstates, sys.noutputs, sys.ninputs, strictly_proper=True)
+        if fcn not in [ct.initial_response, ct.forced_response,
+                       ct.input_output_response] and \
+           isinstance(sys, ct.lti.LTI):
+            # Reuse the previously computed time to make plots look nicer
+            fcn(newsys, *args, T=response.time[-1]).plot(**kwargs)
+        else:
+            # Compute and plot new response (time is one of the arguments)
+            fcn(newsys, *args).plot(**kwargs)
+
+        # Make sure we have the same axes
+        new_axes = plt.gcf().get_axes()
+        assert new_axes == old_axes
+
+        # Make sure every axes has more than one line
+        for ax in new_axes:
+            assert len(ax.get_lines()) > 1
+
+    # Update the title so we can see what is going on
+    fig = out[0, 0][0].axes.figure
+    fig.suptitle(
+        fig._suptitle._text +
+        f" [{sys.noutputs}x{sys.ninputs}, cs={cmbsig}, "
+        f"ct={cmbtrc}, pi={pltinp}, tr={trpose}]",
+        fontsize='small')
+
+    # Get rid of the figure to free up memory
+    if clear:
+        plt.clf()
+
+
+def test_axes_setup():
+    get_plot_axes = ct.timeplot.get_plot_axes
+
+    sys_2x3 = ct.rss(4, 2, 3)
+    sys_2x3b = ct.rss(4, 2, 3)
+    sys_3x2 = ct.rss(4, 3, 2)
+    sys_3x1 = ct.rss(4, 3, 1)
+
+    # Two plots of the same size leaves axes unchanged
+    out1 = ct.step_response(sys_2x3).plot()
+    out2 = ct.step_response(sys_2x3b).plot()
+    np.testing.assert_equal(get_plot_axes(out1), get_plot_axes(out2))
+    plt.close()
+
+    # Two plots of same net size leaves axes unchanged (unfortunately)
+    out1 = ct.step_response(sys_2x3).plot()
+    out2 = ct.step_response(sys_3x2).plot()
+    np.testing.assert_equal(
+        get_plot_axes(out1).reshape(-1), get_plot_axes(out2).reshape(-1))
+    plt.close()
+
+    # Plots of different shapes generate new plots
+    out1 = ct.step_response(sys_2x3).plot()
+    out2 = ct.step_response(sys_3x1).plot()
+    ax1_list = get_plot_axes(out1).reshape(-1).tolist()
+    ax2_list = get_plot_axes(out2).reshape(-1).tolist()
+    for ax in ax1_list:
+        assert ax not in ax2_list
+    plt.close()
+
+    # Passing a list of axes preserves those axes
+    out1 = ct.step_response(sys_2x3).plot()
+    out2 = ct.step_response(sys_3x1).plot()
+    out3 = ct.step_response(sys_2x3b).plot(ax=get_plot_axes(out1))
+    np.testing.assert_equal(get_plot_axes(out1), get_plot_axes(out3))
+    plt.close()
+
+    # Sending an axes array of the wrong size raises exception
+    with pytest.raises(ValueError, match="not the right shape"):
+        out = ct.step_response(sys_2x3).plot()
+        ct.step_response(sys_3x1).plot(ax=get_plot_axes(out))
+    sys_2x3 = ct.rss(4, 2, 3)
+    sys_2x3b = ct.rss(4, 2, 3)
+    sys_3x2 = ct.rss(4, 3, 2)
+    sys_3x1 = ct.rss(4, 3, 1)
+
+
+@slycotonly
+def test_legend_map():
+    sys_mimo = ct.tf2ss(
+        [[[1], [0.1]], [[0.2], [1]]],
+        [[[1, 0.6, 1], [1, 1, 1]], [[1, 0.4, 1], [1, 2, 1]]], name="MIMO")
+    response = ct.step_response(sys_mimo)
+    response.plot(
+        legend_map=np.array([['center', 'upper right'],
+                             [None, 'center right']]),
+        plot_inputs=True, overlay_signals=True, transpose=True,
+        title='MIMO step response with custom legend placement')
+
+
+def test_combine_time_responses():
+    sys_mimo = ct.rss(4, 2, 2)
+    timepts = np.linspace(0, 10, 100)
+
+    # Combine two response with ntrace = 0
+    U = np.vstack([np.sin(timepts), np.cos(2*timepts)])
+    resp1 = ct.input_output_response(sys_mimo, timepts, U)
+
+    U = np.vstack([np.cos(2*timepts), np.sin(timepts)])
+    resp2 = ct.input_output_response(sys_mimo, timepts, U)
+
+    combresp1 = ct.combine_time_responses([resp1, resp2])
+    assert combresp1.ntraces == 2
+    np.testing.assert_equal(combresp1.y[:, 0, :], resp1.y)
+    np.testing.assert_equal(combresp1.y[:, 1, :], resp2.y)
+
+    # Combine two responses with ntrace != 0
+    resp3 = ct.step_response(sys_mimo, timepts)
+    resp4 = ct.step_response(sys_mimo, timepts)
+    combresp2 = ct.combine_time_responses([resp3, resp4])
+    assert combresp2.ntraces == resp3.ntraces + resp4.ntraces
+    np.testing.assert_equal(combresp2.y[:, 0:2, :], resp3.y)
+    np.testing.assert_equal(combresp2.y[:, 2:4, :], resp4.y)
+
+    # Mixture
+    combresp3 = ct.combine_time_responses([resp1, resp2, resp3])
+    assert combresp3.ntraces == resp3.ntraces + resp4.ntraces
+    np.testing.assert_equal(combresp3.y[:, 0, :], resp1.y)
+    np.testing.assert_equal(combresp3.y[:, 1, :], resp2.y)
+    np.testing.assert_equal(combresp3.y[:, 2:4, :], resp3.y)
+    assert combresp3.trace_types == [None, None] + resp3.trace_types
+    assert combresp3.trace_labels == \
+        [resp1.title, resp2.title] + resp3.trace_labels
+
+    # Rename the traces
+    labels = ["T1", "T2", "T3", "T4"]
+    combresp4 = ct.combine_time_responses(
+        [resp1, resp2, resp3], trace_labels=labels)
+    assert combresp4.trace_labels == labels
+
+    # Automatically generated trace label names and types
+    resp5 = ct.step_response(sys_mimo, timepts)
+    resp5.title = "test"
+    resp5.trace_labels = None
+    resp5.trace_types = None
+    combresp5 = ct.combine_time_responses([resp1, resp5])
+    assert combresp5.trace_labels == [resp1.title] + \
+        ["test, trace 0", "test, trace 1"]
+    assert combresp4.trace_types == [None, None, 'step', 'step']
+
+    with pytest.raises(ValueError, match="must have the same number"):
+        resp = ct.step_response(ct.rss(4, 2, 3), timepts)
+        combresp = ct.combine_time_responses([resp1, resp])
+
+    with pytest.raises(ValueError, match="trace labels does not match"):
+        combresp = ct.combine_time_responses(
+            [resp1, resp2], trace_labels=["T1", "T2", "T3"])
+
+    with pytest.raises(ValueError, match="must have the same time"):
+        resp = ct.step_response(ct.rss(4, 2, 3), timepts/2)
+        combresp6 = ct.combine_time_responses([resp1, resp])
+
+
+@slycotonly
+def test_linestyles():
+    # Check to make sure we can change line styles
+    sys_mimo = ct.tf2ss(
+        [[[1], [0.1]], [[0.2], [1]]],
+        [[[1, 0.6, 1], [1, 1, 1]], [[1, 0.4, 1], [1, 2, 1]]], name="MIMO")
+    out = ct.step_response(sys_mimo).plot('k--', plot_inputs=True)
+    for ax in np.nditer(out, flags=["refs_ok"]):
+        for line in ax.item():
+            assert line.get_color() == 'k'
+            assert line.get_linestyle() == '--'
+
+    out = ct.step_response(sys_mimo).plot(
+        plot_inputs='overlay', overlay_signals=True, overlay_traces=True,
+        output_props=[{'color': c} for c in ['blue', 'orange']],
+        input_props=[{'color': c} for c in ['red', 'green']],
+        trace_props=[{'linestyle': s} for s in ['-', '--']])
+
+    assert out.shape == (1, 1)
+    lines = out[0, 0]
+    assert lines[0].get_color() == 'blue' and lines[0].get_linestyle() == '-'
+    assert lines[1].get_color() == 'orange' and lines[1].get_linestyle() == '-'
+    assert lines[2].get_color() == 'red' and lines[2].get_linestyle() == '-'
+    assert lines[3].get_color() == 'green' and lines[3].get_linestyle() == '-'
+    assert lines[4].get_color() == 'blue' and lines[4].get_linestyle() == '--'
+    assert lines[5].get_color() == 'orange' and lines[5].get_linestyle() == '--'
+    assert lines[6].get_color() == 'red' and lines[6].get_linestyle() == '--'
+    assert lines[7].get_color() == 'green' and lines[7].get_linestyle() == '--'
+
+
+def test_rcParams():
+    sys = ct.rss(2, 2, 2)
+
+    # Create new set of rcParams
+    my_rcParams = {
+        'axes.labelsize': 10,
+        'axes.titlesize': 10,
+        'figure.titlesize': 12,
+        'legend.fontsize': 10,
+        'xtick.labelsize': 10,
+        'ytick.labelsize': 10,
+    }
+
+    # Generate a figure with the new rcParams
+    out = ct.step_response(sys).plot(rcParams=my_rcParams)
+    ax = out[0, 0][0].axes
+    fig = ax.figure
+
+    # Check to make sure new settings were used
+    assert ax.xaxis.get_label().get_fontsize() == 10
+    assert ax.yaxis.get_label().get_fontsize() == 10
+    assert ax.title.get_fontsize() == 10
+    assert ax.xaxis._get_tick_label_size('x') == 10
+    assert ax.yaxis._get_tick_label_size('y') == 10
+    assert fig._suptitle.get_fontsize() == 12
+
+def test_relabel():
+    sys1 = ct.rss(2, inputs='u', outputs='y')
+    sys2 = ct.rss(1, 1, 1)      # uses default i/o labels
+
+    # Generate a plot with specific labels
+    ct.step_response(sys1).plot()
+
+    # Generate a new plot, which overwrites labels
+    out = ct.step_response(sys2).plot()
+    ax = ct.get_plot_axes(out)
+    assert ax[0, 0].get_ylabel() == 'y[0]'
+
+    # Regenerate the first plot
+    plt.figure()
+    ct.step_response(sys1).plot()
+
+    # Generate a new plt, without relabeling
+    out = ct.step_response(sys2).plot(relabel=False)
+    ax = ct.get_plot_axes(out)
+    assert ax[0, 0].get_ylabel() == 'y'
+
+
+def test_errors():
+    sys = ct.rss(2, 1, 1)
+    stepresp = ct.step_response(sys)
+    with pytest.raises(AttributeError,
+                       match="(has no property|unexpected keyword)"):
+        stepresp.plot(unknown=None)
+
+    with pytest.raises(AttributeError,
+                       match="(has no property|unexpected keyword)"):
+        ct.time_response_plot(stepresp, unknown=None)
+
+    with pytest.raises(ValueError, match="unrecognized value"):
+        stepresp.plot(plot_inputs='unknown')
+
+    for kw in ['input_props', 'output_props', 'trace_props']:
+        propkw = {kw: {'color': 'green'}}
+        with pytest.warns(UserWarning, match="ignored since fmt string"):
+            out = stepresp.plot('k-', **propkw)
+            assert out[0, 0][0].get_color() == 'k'
+
+if __name__ == "__main__":
+    #
+    # Interactive mode: generate plots for manual viewing
+    #
+    # Running this script in python (or better ipython) will show a
+    # collection of figures that should all look OK on the screeen.
+    #
+
+    # In interactive mode, turn on ipython interactive graphics
+    plt.ion()
+
+    # Start by clearing existing figures
+    plt.close('all')
+
+    # Define a set of systems to test
+    sys_siso = ct.tf2ss([1], [1, 2, 1], name="SISO")
+    sys_mimo = ct.tf2ss(
+        [[[1], [0.1]], [[0.2], [1]]],
+        [[[1, 0.6, 1], [1, 1, 1]], [[1, 0.4, 1], [1, 2, 1]]], name="MIMO")
+
+    # Define and run a selected set of interesting tests
+    # def test_response_plots(
+    #      fcn, sys, plot_inputs, plot_outputs, overlay_signals,
+    #      overlay_traces, transpose, second_system, clear=True):
+    N, T, F = None, True, False
+    test_cases = [
+        # response fcn       system    in         out cs ct tr ss
+        (ct.step_response,   sys_siso, N,         T,  F, F, F, F), # 1
+        (ct.step_response,   sys_siso, T,         F,  F, F, F, F), # 2
+        (ct.step_response,   sys_siso, T,         T,  F, F, F, T), # 3
+        (ct.step_response,   sys_siso, 'overlay', T,  F, F, F, T), # 4
+        (ct.step_response,   sys_mimo, F,         T,  F, F, F, F), # 5
+        (ct.step_response,   sys_mimo, T,         T,  F, F, F, F), # 6
+        (ct.step_response,   sys_mimo, 'overlay', T,  F, F, F, F), # 7
+        (ct.step_response,   sys_mimo, T,         T,  T, F, F, F), # 8
+        (ct.step_response,   sys_mimo, T,         T,  T, T, F, F), # 9
+        (ct.step_response,   sys_mimo, T,         T,  F, F, T, F), # 10
+        (ct.step_response,   sys_mimo, T,         T,  T, F, T, F), # 11
+        (ct.step_response,   sys_mimo, 'overlay', T,  T, F, T, F), # 12
+        (ct.forced_response, sys_mimo, N,         T,  T, F, T, F), # 13
+        (ct.forced_response, sys_mimo, 'overlay', T,  F, F, F, F), # 14
+    ]
+    for args in test_cases:
+        test_response_plots(*args, clear=F)
+
+    #
+    # Run a few more special cases to show off capabilities (and save some
+    # of them for use in the documentation).
+    #
+
+    test_legend_map()           # show ability to set legend location
+
+    # Basic step response
+    plt.figure()
+    ct.step_response(sys_mimo).plot()
+    plt.savefig('timeplot-mimo_step-default.png')
+
+    # Step response with plot_inputs, overlay_signals
+    plt.figure()
+    ct.step_response(sys_mimo).plot(
+        plot_inputs=True, overlay_signals=True,
+        title="Step response for 2x2 MIMO system " +
+        "[plot_inputs, overlay_signals]")
+    plt.savefig('timeplot-mimo_step-pi_cs.png')
+
+    # Input/output response with overlaid inputs, legend_map
+    plt.figure()
+    timepts = np.linspace(0, 10, 100)
+    U = np.vstack([np.sin(timepts), np.cos(2*timepts)])
+    ct.input_output_response(sys_mimo, timepts, U).plot(
+        plot_inputs='overlay',
+        legend_map=np.array([['lower right'], ['lower right']]),
+        title="I/O response for 2x2 MIMO system " +
+        "[plot_inputs='overlay', legend_map]")
+    plt.savefig('timeplot-mimo_ioresp-ov_lm.png')
+
+    # Multi-trace plot, transpose
+    plt.figure()
+    U = np.vstack([np.sin(timepts), np.cos(2*timepts)])
+    resp1 = ct.input_output_response(sys_mimo, timepts, U)
+
+    U = np.vstack([np.cos(2*timepts), np.sin(timepts)])
+    resp2 = ct.input_output_response(sys_mimo, timepts, U)
+
+    ct.combine_time_responses(
+        [resp1, resp2], trace_labels=["Scenario #1", "Scenario #2"]).plot(
+            transpose=True,
+            title="I/O responses for 2x2 MIMO system, multiple traces "
+            "[transpose]")
+    plt.savefig('timeplot-mimo_ioresp-mt_tr.png')
+
+    plt.figure()
+    out = ct.step_response(sys_mimo).plot(
+    plot_inputs='overlay', overlay_signals=True, overlay_traces=True,
+    output_props=[{'color': c} for c in ['blue', 'orange']],
+    input_props=[{'color': c} for c in ['red', 'green']],
+    trace_props=[{'linestyle': s} for s in ['-', '--']])
+    plt.savefig('timeplot-mimo_step-linestyle.png')
diff --git a/control/timeplot.py b/control/timeplot.py
new file mode 100644
index 000000000..6409a6660
--- /dev/null
+++ b/control/timeplot.py
@@ -0,0 +1,808 @@
+# timeplot.py - time plotting functions
+# RMM, 20 Jun 2023
+#
+# This file contains routines for plotting out time responses.  These
+# functions can be called either as standalone functions or access from the
+# TimeDataResponse class.
+#
+# Note: It might eventually make sense to put the functions here
+# directly into timeresp.py.
+
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from os.path import commonprefix
+from warnings import warn
+
+from . import config
+
+__all__ = ['time_response_plot', 'combine_time_responses', 'get_plot_axes']
+
+# Default font dictionary
+_timeplot_rcParams = mpl.rcParams.copy()
+_timeplot_rcParams.update({
+    'axes.labelsize': 'small',
+    'axes.titlesize': 'small',
+    'figure.titlesize': 'medium',
+    'legend.fontsize': 'x-small',
+    'xtick.labelsize': 'small',
+    'ytick.labelsize': 'small',
+})
+
+# Default values for module parameter variables
+_timeplot_defaults = {
+    'timeplot.rcParams': _timeplot_rcParams,
+    'timeplot.trace_props': [
+        {'linestyle': s} for s in ['-', '--', ':', '-.']],
+    'timeplot.output_props': [
+        {'color': c} for c in [
+            'tab:blue', 'tab:orange', 'tab:green', 'tab:pink', 'tab:gray']],
+    'timeplot.input_props': [
+        {'color': c} for c in [
+            'tab:red', 'tab:purple', 'tab:brown', 'tab:olive', 'tab:cyan']],
+    'timeplot.time_label': "Time [s]",
+}
+
+# Plot the input/output response of a system
+def time_response_plot(
+        data, *fmt, ax=None, plot_inputs=None, plot_outputs=True,
+        transpose=False, overlay_traces=False, overlay_signals=False,
+        legend_map=None, legend_loc=None, add_initial_zero=True,
+        trace_labels=None, title=None, relabel=True, **kwargs):
+    """Plot the time response of an input/output system.
+
+    This function creates a standard set of plots for the input/output
+    response of a system, with the data provided via a `TimeResponseData`
+    object, which is the standard output for python-control simulation
+    functions.
+
+    Parameters
+    ----------
+    data : TimeResponseData
+        Data to be plotted.
+    ax : array of Axes
+        The matplotlib Axes to draw the figure on.  If not specified, the
+        Axes for the current figure are used or, if there is no current
+        figure with the correct number and shape of Axes, a new figure is
+        created.  The default shape of the array should be (noutputs +
+        ninputs, ntraces), but if `overlay_traces` is set to `True` then
+        only one row is needed and if `overlay_signals` is set to `True`
+        then only one or two columns are needed (depending on plot_inputs
+        and plot_outputs).
+    plot_inputs : bool or str, optional
+        Sets how and where to plot the inputs:
+            * False: don't plot the inputs
+            * None: use value from time response data (default)
+            * 'overlay`: plot inputs overlaid with outputs
+            * True: plot the inputs on their own axes
+    plot_outputs : bool, optional
+        If False, suppress plotting of the outputs.
+    overlay_traces : bool, optional
+        If set to True, combine all traces onto a single row instead of
+        plotting a separate row for each trace.
+    overlay_signals : bool, optional
+        If set to True, combine all input and output signals onto a single
+        plot (for each).
+    transpose : bool, optional
+        If transpose is False (default), signals are plotted from top to
+        bottom, starting with outputs (if plotted) and then inputs.
+        Multi-trace plots are stacked horizontally.  If transpose is True,
+        signals are plotted from left to right, starting with the inputs
+        (if plotted) and then the outputs.  Multi-trace responses are
+        stacked vertically.
+    *fmt : :func:`matplotlib.pyplot.plot` format string, optional
+        Passed to `matplotlib` as the format string for all lines in the plot.
+    **kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional
+        Additional keywords passed to `matplotlib` to specify line properties.
+
+    Returns
+    -------
+    out : array of list of Line2D
+        Array of Line2D objects for each line in the plot.  The shape of
+        the array matches the subplots shape and the value of the array is a
+        list of Line2D objects in that subplot.
+
+    Additional Parameters
+    ---------------------
+    add_initial_zero : bool
+        Add an initial point of zero at the first time point for all
+        inputs with type 'step'.  Default is True.
+    input_props : array of dicts
+        List of line properties to use when plotting combined inputs.  The
+        default values are set by config.defaults['timeplot.input_props'].
+    legend_map : array of str, option
+        Location of the legend for multi-trace plots.  Specifies an array
+        of legend location strings matching the shape of the subplots, with
+        each entry being either None (for no legend) or a legend location
+        string (see :func:`~matplotlib.pyplot.legend`).
+    legend_loc : str
+        Location of the legend within the axes for which it appears.  This
+        value is used if legend_map is None.
+    output_props : array of dicts
+        List of line properties to use when plotting combined outputs.  The
+        default values are set by config.defaults['timeplot.output_props'].
+    relabel : bool, optional
+        By default, existing figures and axes are relabeled when new data
+        are added.  If set to `False`, just plot new data on existing axes.
+    time_label : str, optional
+        Label to use for the time axis.
+    trace_props : array of dicts
+        List of line properties to use when plotting combined outputs.  The
+        default values are set by config.defaults['timeplot.trace_props'].
+
+    Notes
+    -----
+    1. A new figure will be generated if there is no current figure or
+       the current figure has an incompatible number of axes.  To
+       force the creation of a new figures, use `plt.figure()`.  To reuse
+       a portion of an existing figure, use the `ax` keyword.
+
+    2. The line properties (color, linestyle, etc) can be set for the
+       entire plot using the `fmt` and/or `kwargs` parameter, which
+       are passed on to `matplotlib`.  When combining signals or
+       traces, the `input_props`, `output_props`, and `trace_props`
+       parameters can be used to pass a list of dictionaries
+       containing the line properties to use.  These input/output
+       properties are combined with the trace properties and finally
+       the kwarg properties to determine the final line properties.
+
+    3. The default plot properties, such as font sizes, can be set using
+       config.defaults[''timeplot.rcParams'].
+
+    """
+    from .iosys import InputOutputSystem
+    from .timeresp import TimeResponseData
+
+    #
+    # Process keywords and set defaults
+    #
+
+    # Set up defaults
+    time_label = config._get_param(
+        'timeplot', 'time_label', kwargs, _timeplot_defaults, pop=True)
+    timeplot_rcParams = config._get_param(
+        'timeplot', 'rcParams', kwargs, _timeplot_defaults, pop=True)
+
+    if kwargs.get('input_props', None) and len(fmt) > 0:
+        warn("input_props ignored since fmt string was present")
+    input_props = config._get_param(
+        'timeplot', 'input_props', kwargs, _timeplot_defaults, pop=True)
+    iprop_len = len(input_props)
+
+    if kwargs.get('output_props', None) and len(fmt) > 0:
+        warn("output_props ignored since fmt string was present")
+    output_props = config._get_param(
+        'timeplot', 'output_props', kwargs, _timeplot_defaults, pop=True)
+    oprop_len = len(output_props)
+
+    if kwargs.get('trace_props', None) and len(fmt) > 0:
+        warn("trace_props ignored since fmt string was present")
+    trace_props = config._get_param(
+        'timeplot', 'trace_props', kwargs, _timeplot_defaults, pop=True)
+    tprop_len = len(trace_props)
+
+    # Set the title for the data
+    title = data.title if title == None else title
+
+    # Determine whether or not to plot the input data (and how)
+    if plot_inputs is None:
+        plot_inputs = data.plot_inputs
+    if plot_inputs not in [True, False, 'overlay']:
+        raise ValueError(f"unrecognized value: {plot_inputs=}")
+
+    #
+    # Find/create axes
+    #
+    # Data are plotted in a standard subplots array, whose size depends on
+    # which signals are being plotted and how they are combined.  The
+    # baseline layout for data is to plot everything separately, with
+    # outputs and inputs making up the rows and traces making up the
+    # columns:
+    #
+    # Trace 0        Trace q
+    # +------+       +------+
+    # | y[0] |  ...  | y[0] |
+    # +------+       +------+
+    #    :
+    # +------+       +------+
+    # | y[p] |  ...  | y[p] |
+    # +------+       +------+
+    #
+    # +------+       +------+
+    # | u[0] |  ...  | u[0] |
+    # +------+       +------+
+    #    :
+    # +------+       +------+
+    # | u[m] |  ...  | u[m] |
+    # +------+       +------+
+    #
+    # A variety of options are available to modify this format:
+    #
+    # * Omitting: either the inputs or the outputs can be omitted.
+    #
+    # * Combining: inputs, outputs, and traces can be combined onto a
+    #   single set of axes using various keyword combinations
+    #   (overlay_signals, overlay_traces, plot_inputs='overlay').  This
+    #   basically collapses data along either the rows or columns, and a
+    #   legend is generated.
+    #
+    # * Transpose: if the `transpose` keyword is True, then instead of
+    #   plotting the data vertically (outputs over inputs), we plot left to
+    #   right (inputs, outputs):
+    #
+    #         +------+       +------+     +------+       +------+
+    # Trace 0 | u[0] |  ...  | u[m] |     | y[0] |  ...  | y[p] |
+    #         +------+       +------+     +------+       +------+
+    #    :
+    #    :
+    #         +------+       +------+     +------+       +------+
+    # Trace q | u[0] |  ...  | u[m] |     | y[0] |  ...  | y[p] |
+    #         +------+       +------+     +------+       +------+
+    #
+    # This also affects the way in which legends and labels are generated.
+
+    # Decide on the number of inputs and outputs
+    ninputs = data.ninputs if plot_inputs else 0
+    noutputs = data.noutputs if plot_outputs else 0
+    ntraces = max(1, data.ntraces)      # treat data.ntraces == 0 as 1 trace
+    if ninputs == 0 and noutputs == 0:
+        raise ValueError(
+            "plot_inputs and plot_outputs both False; no data to plot")
+    elif plot_inputs == 'overlay' and noutputs == 0:
+        raise ValueError(
+            "can't overlay inputs with no outputs")
+    elif plot_inputs in [True, 'overlay'] and data.ninputs == 0:
+        raise ValueError(
+            "input plotting requested but no inputs in time response data")
+
+    # Figure how how many rows and columns to use + offsets for inputs/outputs
+    if plot_inputs == 'overlay' and not overlay_signals:
+        nrows = max(ninputs, noutputs)          # Plot inputs on top of outputs
+        noutput_axes = 0                        # No offset required
+        ninput_axes = 0                         # No offset required
+    elif overlay_signals:
+        nrows = int(plot_outputs)               # Start with outputs
+        nrows += int(plot_inputs == True)       # Add plot for inputs if needed
+        noutput_axes = 1 if plot_outputs and plot_inputs is True else 0
+        ninput_axes = 1 if plot_inputs is True else 0
+    else:
+        nrows = noutputs + ninputs              # Plot inputs separately
+        noutput_axes = noutputs if plot_outputs else 0
+        ninput_axes = ninputs if plot_inputs else 0
+
+    ncols = ntraces if not overlay_traces else 1
+    if transpose:
+        nrows, ncols = ncols, nrows
+
+    # See if we can use the current figure axes
+    fig = plt.gcf()         # get current figure (or create new one)
+    if ax is None and plt.get_fignums():
+        ax = fig.get_axes()
+        if len(ax) == nrows * ncols:
+            # Assume that the shape is right (no easy way to infer this)
+            ax = np.array(ax).reshape(nrows, ncols)
+        elif len(ax) != 0:
+            # Need to generate a new figure
+            fig, ax = plt.figure(), None
+        else:
+            # Blank figure, just need to recreate axes
+            ax = None
+
+    # Create new axes, if needed, and customize them
+    if ax is None:
+        with plt.rc_context(timeplot_rcParams):
+            ax_array = fig.subplots(nrows, ncols, sharex=True, squeeze=False)
+            fig.set_tight_layout(True)
+            fig.align_labels()
+
+    else:
+        # Make sure the axes are the right shape
+        if ax.shape != (nrows, ncols):
+            raise ValueError(
+                "specified axes are not the right shape; "
+                f"got {ax.shape} but expecting ({nrows}, {ncols})")
+        ax_array = ax
+
+    #
+    # Map inputs/outputs and traces to axes
+    #
+    # This set of code takes care of all of the various options for how to
+    # plot the data.  The arrays output_map and input_map are used to map
+    # the different signals that are plotted onto the axes created above.
+    # This code is complicated because it has to handle lots of different
+    # variations.
+    #
+
+    # Create the map from trace, signal to axes, accounting for overlay_*
+    output_map = np.empty((noutputs, ntraces), dtype=tuple)
+    input_map = np.empty((ninputs, ntraces), dtype=tuple)
+
+    for i in range(noutputs):
+        for j in range(ntraces):
+            signal_index = i if not overlay_signals else 0
+            trace_index = j if not overlay_traces else 0
+            if transpose:
+                output_map[i, j] = (trace_index, signal_index + ninput_axes)
+            else:
+                output_map[i, j] = (signal_index, trace_index)
+
+    for i in range(ninputs):
+        for j in range(ntraces):
+            signal_index = noutput_axes + (i if not overlay_signals else 0)
+            trace_index = j if not overlay_traces else 0
+            if transpose:
+                input_map[i, j] = (trace_index, signal_index - noutput_axes)
+            else:
+                input_map[i, j] = (signal_index, trace_index)
+
+    #
+    # Plot the data
+    #
+    # The ax_output and ax_input arrays have the axes needed for making the
+    # plots.  Labels are used on each axes for later creation of legends.
+    # The gneric labels if of the form:
+    #
+    #     signal name, trace label, system name
+    #
+    # The signal name or tracel label can be omitted if they will appear on
+    # the axes title or ylabel.  The system name is always included, since
+    # multiple calls to plot() will require a legend that distinguishes
+    # which system signals are plotted.  The system name is stripped off
+    # later (in the legend-handling code) if it is not needed, but must be
+    # included here since a plot may be built up by multiple calls to plot().
+    #
+
+    # Reshape the inputs and outputs for uniform indexing
+    outputs = data.y.reshape(data.noutputs, ntraces, -1)
+    if data.u is None or not plot_inputs:
+        inputs = None
+    else:
+        inputs = data.u.reshape(data.ninputs, ntraces, -1)
+
+    # Create a list of lines for the output
+    out = np.empty((nrows, ncols), dtype=object)
+    for i in range(nrows):
+        for j in range(ncols):
+            out[i, j] = []      # unique list in each element
+
+    # Utility function for creating line label
+    def _make_line_label(signal_index, signal_labels, trace_index):
+        label = ""              # start with an empty label
+
+        # Add the signal name if it won't appear as an axes label
+        if overlay_signals or plot_inputs == 'overlay':
+            label += signal_labels[signal_index]
+
+        # Add the trace label if this is a multi-trace figure
+        if overlay_traces and ntraces > 1 or trace_labels:
+            label += ", " if label != "" else ""
+            if trace_labels:
+                label += trace_labels[trace_index]
+            elif data.trace_labels:
+                label += data.trace_labels[trace_index]
+            else:
+                label += f"trace {trace_index}"
+
+        # Add the system name (will strip off later if redundant)
+        label += ", " if label != "" else ""
+        label += f"{data.sysname}"
+
+        return label
+
+    # Go through each trace and each input/output
+    for trace in range(ntraces):
+        # Plot the output
+        for i in range(noutputs):
+            label = _make_line_label(i, data.output_labels, trace)
+
+            # Set up line properties for this output, trace
+            if len(fmt) == 0:
+                line_props = output_props[
+                    i % oprop_len if overlay_signals else 0].copy()
+                line_props.update(
+                    trace_props[trace % tprop_len if overlay_traces else 0])
+                line_props.update(kwargs)
+            else:
+                line_props = kwargs
+
+            out[output_map[i, trace]] += ax_array[output_map[i, trace]].plot(
+                data.time, outputs[i][trace], *fmt, label=label, **line_props)
+
+        # Plot the input
+        for i in range(ninputs):
+            label = _make_line_label(i, data.input_labels, trace)
+
+            if add_initial_zero and data.ntraces > i \
+               and data.trace_types[i] == 'step':
+                x = np.hstack([np.array([data.time[0]]), data.time])
+                y = np.hstack([np.array([0]), inputs[i][trace]])
+            else:
+                x, y = data.time, inputs[i][trace]
+
+            # Set up line properties for this output, trace
+            if len(fmt) == 0:
+                line_props = input_props[
+                    i % iprop_len if overlay_signals else 0].copy()
+                line_props.update(
+                    trace_props[trace % tprop_len if overlay_traces else 0])
+                line_props.update(kwargs)
+            else:
+                line_props = kwargs
+
+            out[input_map[i, trace]] += ax_array[input_map[i, trace]].plot(
+                x, y, *fmt, label=label, **line_props)
+
+    # Stop here if the user wants to control everything
+    if not relabel:
+        return out
+
+    #
+    # Label the axes (including trace labels)
+    #
+    # Once the data are plotted, we label the axes.  The horizontal axes is
+    # always time and this is labeled only on the bottom most column.  The
+    # vertical axes can consist either of a single signal or a combination
+    # of signals (when overlay_signal is True or plot+inputs = 'overlay'.
+    #
+    # Traces are labeled at the top of the first row of plots (regular) or
+    # the left edge of rows (tranpose).
+    #
+
+    # Time units on the bottom
+    for col in range(ncols):
+        ax_array[-1, col].set_xlabel(time_label)
+
+    # Keep track of whether inputs are overlaid on outputs
+    overlaid = plot_inputs == 'overlay'
+    overlaid_title = "Inputs, Outputs"
+
+    if transpose:               # inputs on left, outputs on right
+        # Label the inputs
+        if overlay_signals and plot_inputs:
+            label = overlaid_title if overlaid else "Inputs"
+            for trace in range(ntraces):
+                ax_array[input_map[0, trace]].set_ylabel(label)
+        else:
+            for i in range(ninputs):
+                label = overlaid_title if overlaid else data.input_labels[i]
+                for trace in range(ntraces):
+                    ax_array[input_map[i, trace]].set_ylabel(label)
+
+        # Label the outputs
+        if overlay_signals and plot_outputs:
+            label = overlaid_title if overlaid else "Outputs"
+            for trace in range(ntraces):
+                ax_array[output_map[0, trace]].set_ylabel(label)
+        else:
+            for i in range(noutputs):
+                label = overlaid_title if overlaid else data.output_labels[i]
+                for trace in range(ntraces):
+                    ax_array[output_map[i, trace]].set_ylabel(label)
+
+        # Set the trace titles, if needed
+        if ntraces > 1 and not overlay_traces:
+            for trace in range(ntraces):
+                # Get the existing ylabel for left column
+                label = ax_array[trace, 0].get_ylabel()
+
+                # Add on the trace title
+                if trace_labels:
+                    label = trace_labels[trace] + "\n" + label
+                elif data.trace_labels:
+                    label = data.trace_labels[trace] + "\n" + label
+                else:
+                    label = f"Trace {trace}" + "\n" + label
+
+                ax_array[trace, 0].set_ylabel(label)
+
+    else:                       # regular plot (outputs over inputs)
+        # Set the trace titles, if needed
+        if ntraces > 1 and not overlay_traces:
+            for trace in range(ntraces):
+                if trace_labels:
+                    label = trace_labels[trace]
+                elif data.trace_labels:
+                    label = data.trace_labels[trace]
+                else:
+                    label = f"Trace {trace}"
+
+                with plt.rc_context(timeplot_rcParams):
+                    ax_array[0, trace].set_title(label)
+
+        # Label the outputs
+        if overlay_signals and plot_outputs:
+            ax_array[output_map[0, 0]].set_ylabel("Outputs")
+        else:
+            for i in range(noutputs):
+                ax_array[output_map[i, 0]].set_ylabel(
+                    overlaid_title if overlaid else data.output_labels[i])
+
+        # Label the inputs
+        if overlay_signals and plot_inputs:
+            label = overlaid_title if overlaid else "Inputs"
+            ax_array[input_map[0, 0]].set_ylabel(label)
+        else:
+            for i in range(ninputs):
+                label = overlaid_title if overlaid else data.input_labels[i]
+                ax_array[input_map[i, 0]].set_ylabel(label)
+
+    #
+    # Create legends
+    #
+    # Legends can be placed manually by passing a legend_map array that
+    # matches the shape of the suplots, with each item being a string
+    # indicating the location of the legend for that axes (or None for no
+    # legend).
+    #
+    # If no legend spec is passed, a minimal number of legends are used so
+    # that each line in each axis can be uniquely identified.  The details
+    # depends on the various plotting parameters, but the general rule is
+    # to place legends in the top row and right column.
+    #
+    # Because plots can be built up by multiple calls to plot(), the legend
+    # strings are created from the line labels manually.  Thus an initial
+    # call to plot() may not generate any legends (eg, if no signals are
+    # combined nor overlaid), but subsequent calls to plot() will need a
+    # legend for each different line (system).
+    #
+
+    # Figure out where to put legends
+    if legend_map is None:
+        legend_map = np.full(ax_array.shape, None, dtype=object)
+        if legend_loc == None:
+            legend_loc = 'center right'
+        if transpose:
+            if (overlay_signals or plot_inputs == 'overlay') and overlay_traces:
+                # Put a legend in each plot for inputs and outputs
+                if plot_outputs is True:
+                    legend_map[0, ninput_axes] = legend_loc
+                if plot_inputs is True:
+                    legend_map[0, 0] = legend_loc
+            elif overlay_signals:
+                # Put a legend in rightmost input/output plot
+                if plot_inputs is True:
+                    legend_map[0, 0] = legend_loc
+                if plot_outputs is True:
+                    legend_map[0, ninput_axes] = legend_loc
+            elif plot_inputs == 'overlay':
+                # Put a legend on the top of each column
+                for i in range(ntraces):
+                    legend_map[0, i] = legend_loc
+            elif overlay_traces:
+                # Put a legend topmost input/output plot
+                legend_map[0, -1] = legend_loc
+            else:
+                # Put legend in the upper right
+                legend_map[0, -1] = legend_loc
+        else:                   # regular layout
+            if (overlay_signals or plot_inputs == 'overlay') and overlay_traces:
+                # Put a legend in each plot for inputs and outputs
+                if plot_outputs is True:
+                    legend_map[0, -1] = legend_loc
+                if plot_inputs is True:
+                    legend_map[noutput_axes, -1] = legend_loc
+            elif overlay_signals:
+                # Put a legend in rightmost input/output plot
+                if plot_outputs is True:
+                    legend_map[0, -1] = legend_loc
+                if plot_inputs is True:
+                    legend_map[noutput_axes, -1] = legend_loc
+            elif plot_inputs == 'overlay':
+                # Put a legend on the right of each row
+                for i in range(max(ninputs, noutputs)):
+                    legend_map[i, -1] = legend_loc
+            elif overlay_traces:
+                # Put a legend topmost input/output plot
+                legend_map[0, -1] = legend_loc
+            else:
+                # Put legend in the upper right
+                legend_map[0, -1] = legend_loc
+
+    # Create axis legends
+    for i in range(nrows):
+        for j in range(ncols):
+            ax = ax_array[i, j]
+            # Get the labels to use
+            labels = [line.get_label() for line in ax.get_lines()]
+
+            # Look for a common prefix (up to a space)
+            common_prefix = commonprefix(labels)
+            last_space = common_prefix.rfind(', ')
+            if last_space < 0 or plot_inputs == 'overlay':
+                common_prefix = ''
+            elif last_space > 0:
+                common_prefix = common_prefix[:last_space]
+            prefix_len = len(common_prefix)
+
+            # Look for a common suffice (up to a space)
+            common_suffix = commonprefix(
+                [label[::-1] for label in labels])[::-1]
+            suffix_len = len(common_suffix)
+            # Only chop things off after a comma or space
+            while suffix_len > 0 and common_suffix[-suffix_len] != ',':
+                suffix_len -= 1
+
+            # Strip the labels of common information
+            if suffix_len > 0:
+                labels = [label[prefix_len:-suffix_len] for label in labels]
+            else:
+                labels = [label[prefix_len:] for label in labels]
+
+            # Update the labels to remove common strings
+            if len(labels) > 1 and legend_map[i, j] != None:
+                with plt.rc_context(timeplot_rcParams):
+                    ax.legend(labels, loc=legend_map[i, j])
+
+    #
+    # Update the plot title (= figure suptitle)
+    #
+    # If plots are built up by multiple calls to plot() and the title is
+    # not given, then the title is updated to provide a list of unique text
+    # items in each successive title.  For data generated by the I/O
+    # response functions this will generate a common prefix followed by a
+    # list of systems (e.g., "Step response for sys[1], sys[2]").
+    #
+
+    if fig is not None and title is not None:
+        # Get the current title, if it exists
+        old_title = None if fig._suptitle is None else fig._suptitle._text
+        new_title = title
+
+        if old_title is not None:
+            # Find the common part of the titles
+            common_prefix = commonprefix([old_title, new_title])
+
+            # Back up to the last space
+            last_space = common_prefix.rfind(' ')
+            if last_space > 0:
+                common_prefix = common_prefix[:last_space]
+            common_len = len(common_prefix)
+
+            # Add the new part of the title (usually the system name)
+            if old_title[common_len:] != new_title[common_len:]:
+                separator = ',' if len(common_prefix) > 0 else ';'
+                new_title = old_title + separator + new_title[common_len:]
+
+        # Add the title
+        with plt.rc_context(timeplot_rcParams):
+            fig.suptitle(new_title)
+
+    return out
+
+
+def combine_time_responses(response_list, trace_labels=None, title=None):
+    """Combine multiple individual time responses into a multi-trace response.
+
+    This function combines multiple instances of :class:`TimeResponseData`
+    into a multi-trace :class:`TimeResponseData` object.
+
+    Parameters
+    ----------
+    response_list : list of :class:`TimeResponseData` objects
+        Reponses to be combined.
+    trace_labels : list of str, optional
+        List of labels for each trace.  If not specified, trace names are
+        taken from the input data or set to None.
+
+    Returns
+    -------
+    data : :class:`TimeResponseData`
+        Multi-trace input/output data.
+
+    """
+    from .timeresp import TimeResponseData
+
+    # Save the first trace as the base case
+    base = response_list[0]
+
+    # Process keywords
+    title = base.title if title is None else title
+
+    # Figure out the size of the data (and check for consistency)
+    ntraces = max(1, base.ntraces)
+
+    # Initial pass through trace list to count things up and do error checks
+    for response in response_list[1:]:
+        # Make sure the time vector is the same
+        if not np.allclose(base.t, response.t):
+            raise ValueError("all responses must have the same time vector")
+
+        # Make sure the dimensions are all the same
+        if base.ninputs != response.ninputs or \
+           base.noutputs != response.noutputs or \
+           base.nstates != response.nstates:
+            raise ValueError("all responses must have the same number of "
+                            "inputs, outputs, and states")
+
+        ntraces += max(1, response.ntraces)
+
+    # Create data structures for the new time response data object
+    inputs = np.empty((base.ninputs, ntraces, base.t.size))
+    outputs = np.empty((base.noutputs, ntraces, base.t.size))
+    states = np.empty((base.nstates, ntraces, base.t.size))
+
+    # See whether we should create labels or not
+    if trace_labels is None:
+        generate_trace_labels = True
+        trace_labels = []
+    elif len(trace_labels) != ntraces:
+        raise ValueError(
+            "number of trace labels does not match number of traces")
+    else:
+        generate_trace_labels = False
+
+    offset = 0
+    trace_types = []
+    for response in response_list:
+        if response.ntraces == 0:
+            # Single trace
+            inputs[:, offset, :] = response.u
+            outputs[:, offset, :] = response.y
+            states[:, offset, :] = response.x
+            offset += 1
+
+            # Add on trace label and trace type
+            if generate_trace_labels:
+                trace_labels.append(response.title)
+            trace_types.append(
+                None if response.trace_types is None else response.types[0])
+
+        else:
+            # Save the data
+            for i in range(response.ntraces):
+                inputs[:, offset, :] = response.u[:, i, :]
+                outputs[:, offset, :] = response.y[:, i, :]
+                states[:, offset, :] = response.x[:, i, :]
+
+                # Save the trace labels
+                if generate_trace_labels:
+                    if response.trace_labels is not None:
+                        trace_labels.append(response.trace_labels[i])
+                    else:
+                        trace_labels.append(response.title + f", trace {i}")
+
+                offset += 1
+
+            # Save the trace types
+            if response.trace_types is not None:
+                trace_types += response.trace_types
+            else:
+                trace_types += [None] * response.ntraces
+
+    return TimeResponseData(
+        base.t, outputs, states, inputs, issiso=base.issiso,
+        output_labels=base.output_labels, input_labels=base.input_labels,
+        state_labels=base.state_labels, title=title, transpose=base.transpose,
+        return_x=base.return_x, squeeze=base.squeeze, sysname=base.sysname,
+        trace_labels=trace_labels, trace_types=trace_types,
+        plot_inputs=base.plot_inputs)
+
+
+# Create vectorized function to find axes from lines
+def get_plot_axes(line_array):
+    """Get a list of axes from an array of lines.
+
+    This function can be used to return the set of axes corresponding to
+    the line array that is returned by `time_response_plot`.  This is useful for
+    generating an axes array that can be passed to subsequent plotting
+    calls.
+
+    Parameters
+    ----------
+    line_array : array of list of Line2D
+        A 2D array with elements corresponding to a list of lines appearing
+        in an axes, matching the return type of a time response data plot.
+
+    Returns
+    -------
+    axes_array : array of list of Axes
+        A 2D array with elements corresponding to the Axes assocated with
+        the lines in `line_array`.
+
+    Notes
+    -----
+    Only the first element of each array entry is used to determine the axes.
+
+    """
+    _get_axes = np.vectorize(lambda lines: lines[0].axes)
+    return _get_axes(line_array)
diff --git a/control/timeresp.py b/control/timeresp.py
index defdbdf4e..bc13ad0d1 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -81,6 +81,7 @@
 from . import config
 from .exception import pandas_check
 from .iosys import isctime, isdtime
+from .timeplot import time_response_plot
 
 
 __all__ = ['forced_response', 'step_response', 'step_info',
@@ -169,10 +170,24 @@ class TimeResponseData:
     input_labels, output_labels, state_labels : array of str
         Names for the input, output, and state variables.
 
-    ntraces : int
+    sysname : str, optional
+        Name of the system that created the data.
+
+    plot_inputs : bool, optional
+        Whether or not to plot the inputs by default (can be overridden in
+        the plot() method)
+
+    ntraces : int, optional
         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.
+        response.  If ntraces is 0 (default) then the data represents a
+        single trace with the trace index surpressed in the data.
+
+    trace_labels : array of string, optional
+        Labels to use for traces (set to sysname it ntraces is 0)
+
+    trace_types : array of string, optional
+        Type of trace.  Currently only 'step' is supported, which controls
+        the way in which the signal is plotted.
 
     Notes
     -----
@@ -210,7 +225,9 @@ class TimeResponseData:
     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
+            title=None, transpose=False, return_x=False, squeeze=None,
+            multi_trace=False, trace_labels=None, trace_types=None,
+            plot_inputs=True, sysname=None
     ):
         """Create an input/output time response object.
 
@@ -241,9 +258,8 @@ def __init__(
             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.
+        title : str, optonal
+            Title of the data set (used as figure title in plotting).
 
         squeeze : bool, optional
             By default, if a system is single-input, single-output (SISO)
@@ -267,6 +283,9 @@ def __init__(
             Optional labels for the inputs, outputs, and states, given as a
             list of strings matching the appropriate signal dimension.
 
+        sysname : str, optional
+            Name of the system that created the data.
+
         transpose : bool, optional
             If True, transpose all input and output arrays (for backward
             compatibility with MATLAB and :func:`scipy.signal.lsim`).
@@ -276,6 +295,10 @@ def __init__(
             If True, return the state vector when enumerating result by
             assigning to a tuple (default = False).
 
+        plot_inputs : bool, optional
+            Whether or not to plot the inputs by default (can be overridden
+            in the plot() method)
+
         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
@@ -290,6 +313,8 @@ def __init__(
         self.t = np.atleast_1d(time)
         if self.t.ndim != 1:
             raise ValueError("Time vector must be 1D array")
+        self.title = title
+        self.sysname = sysname
 
         #
         # Output vector (and number of traces)
@@ -363,9 +388,11 @@ def __init__(
         if inputs is None:
             self.u = None
             self.ninputs = 0
+            self.plot_inputs = False
 
         else:
             self.u = np.array(inputs)
+            self.plot_inputs = plot_inputs
 
             # Make sure the shape is OK and figure out the nuumber of inputs
             if multi_trace and self.u.ndim == 3 and \
@@ -397,6 +424,11 @@ def __init__(
         self.input_labels = _process_labels(
             input_labels, "input", self.ninputs)
 
+        # Check and store trace labels, if present
+        self.trace_labels = _process_labels(
+            trace_labels, "trace", self.ntraces)
+        self.trace_types = trace_types
+
         # Figure out if the system is SISO
         if issiso is None:
             # Figure out based on the data
@@ -655,6 +687,10 @@ def to_pandas(self):
 
         return pandas.DataFrame(data)
 
+    # Plot data
+    def plot(self, *args, **kwargs):
+        return time_response_plot(self, *args, **kwargs)
+
 
 # Process signal labels
 def _process_labels(labels, signal, length):
@@ -1132,7 +1168,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     return TimeResponseData(
         tout, yout, xout, U, issiso=sys.issiso(),
         output_labels=sys.output_labels, input_labels=sys.input_labels,
-        state_labels=sys.state_labels,
+        state_labels=sys.state_labels, sysname=sys.name, plot_inputs=True,
+        title="Forced response for " + sys.name, trace_types=['forced'],
         transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
@@ -1324,7 +1361,7 @@ def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
     from .statesp import _convert_to_statespace
 
     # Create the time and input vectors
-    if T is None or np.asarray(T).size == 1 and isinstance(sys, LTI):
+    if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
     T = np.atleast_1d(T).reshape(-1)
     if T.ndim != 1 and len(T) < 2:
@@ -1338,7 +1375,7 @@ def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
             "with given X0.")
 
     # Convert to state space so that we can simulate
-    if sys.nstates is None:
+    if isinstance(sys, LTI) and sys.nstates is None:
         sys = _convert_to_statespace(sys)
 
     # Set up arrays to handle the output
@@ -1349,11 +1386,16 @@ def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
     uout = np.empty((ninputs, ninputs, T.size))
 
     # Simulate the response for each input
+    trace_labels, trace_types = [], []
     for i in range(sys.ninputs):
         # If input keyword was specified, only simulate for that input
         if isinstance(input, int) and i != input:
             continue
 
+        # Save a label and type for this plot
+        trace_labels.append(f"From {sys.input_labels[i]}")
+        trace_types.append('step')
+
         # Create a set of single inputs system for simulation
         U = np.zeros((sys.ninputs, T.size))
         U[i, :] = np.ones_like(T)
@@ -1377,8 +1419,10 @@ def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
     return TimeResponseData(
         response.time, yout, xout, uout, issiso=issiso,
         output_labels=output_labels, input_labels=input_labels,
-        state_labels=sys.state_labels,
-        transpose=transpose, return_x=return_x, squeeze=squeeze)
+        state_labels=sys.state_labels, title="Step response for " + sys.name,
+        transpose=transpose, return_x=return_x, squeeze=squeeze,
+        sysname=sys.name, trace_labels=trace_labels,
+        trace_types=trace_types, plot_inputs=False)
 
 
 def step_info(sysdata, T=None, T_num=None, yfinal=None,
@@ -1695,10 +1739,7 @@ def initial_response(sys, T=None, X0=0, output=None, T_num=None,
 
     # Create the time and input vectors
     if T is None or np.asarray(T).size == 1:
-        if isinstance(sys, LTI):
-            T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)
-        elif T_num is not None:
-            T = np.linspace(0, T, T_num)
+        T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)
     T = np.atleast_1d(T).reshape(-1)
     if T.ndim != 1 and len(T) < 2:
         raise ValueError("invalid value of T for this type of system")
@@ -1718,7 +1759,8 @@ def initial_response(sys, T=None, X0=0, output=None, T_num=None,
     return TimeResponseData(
         response.t, yout, response.x, None, issiso=issiso,
         output_labels=output_labels, input_labels=None,
-        state_labels=sys.state_labels,
+        state_labels=sys.state_labels, sysname=sys.name,
+        title="Initial response for " + sys.name, trace_types=['initial'],
         transpose=transpose, return_x=return_x, squeeze=squeeze)
 
 
@@ -1819,7 +1861,7 @@ def impulse_response(sys, T=None, input=None, output=None, T_num=None,
     from .lti import LTI
 
     # Create the time and input vectors
-    if T is None or np.asarray(T).size == 1 and isinstance(sys, LTI):
+    if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)
     T = np.atleast_1d(T).reshape(-1)
     if T.ndim != 1 and len(T) < 2:
@@ -1845,11 +1887,16 @@ def impulse_response(sys, T=None, input=None, output=None, T_num=None,
     uout = np.full((ninputs, ninputs, np.asarray(T).size), None)
 
     # Simulate the response for each input
+    trace_labels, trace_types = [], []
     for i in range(sys.ninputs):
         # If input keyword was specified, only handle that case
         if isinstance(input, int) and i != input:
             continue
 
+        # Save a label for this plot
+        trace_labels.append(f"From {sys.input_labels[i]}")
+        trace_types.append('impulse')
+
         #
         # Compute new X0 that contains the impulse
         #
@@ -1887,8 +1934,10 @@ def impulse_response(sys, T=None, input=None, output=None, T_num=None,
     return TimeResponseData(
         response.time, yout, xout, uout, issiso=issiso,
         output_labels=output_labels, input_labels=input_labels,
-        state_labels=sys.state_labels,
-        transpose=transpose, return_x=return_x, squeeze=squeeze)
+        state_labels=sys.state_labels, trace_labels=trace_labels,
+        trace_types=trace_types, title="Impulse response for " + sys.name,
+        sysname=sys.name, plot_inputs=False, transpose=transpose,
+        return_x=return_x, squeeze=squeeze)
 
 
 # utility function to find time period and time increment using pole locations
@@ -2062,6 +2111,20 @@ def _ideal_tfinal_and_dt(sys, is_step=True):
 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 """
+    from .lti import LTI
+
+    # For non-LTI system, need tfinal
+    if not isinstance(sys, LTI):
+        if tfinal is None:
+            raise ValueError(
+                "can't automatically compute T for non-LTI system")
+        elif isinstance(tfinal, (int, float, np.number)):
+            if N is None:
+                return np.linspace(0, tfinal)
+            else:
+                return np.linspace(0, tfinal, N)
+        else:
+            return tfinal       # Assume we got passed something appropriate
 
     N_max = 5000
     N_min_ct = 100    # min points for cont time systems
diff --git a/doc/Makefile b/doc/Makefile
index 6e1012343..88a1b7bad 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -15,10 +15,13 @@ help:
 .PHONY: help Makefile
 
 # Rules to create figures
-FIGS = classes.pdf
+FIGS = classes.pdf timeplot-mimo_step-pi_cs.png
 classes.pdf: classes.fig
 	fig2dev -Lpdf $< $@
 
+timeplot-mimo_step-pi_cs.png: ../control/tests/timeplot_test.py
+	PYTHONPATH=.. python $<
+
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
 html pdf clean doctest: Makefile $(FIGS)
diff --git a/doc/index.rst b/doc/index.rst
index 98b184286..ec556e7ce 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -26,6 +26,7 @@ implements basic operations for analysis and design of feedback control systems.
    conventions
    control
    classes
+   plotting
    matlab
    flatsys
    iosys
diff --git a/doc/plotting.rst b/doc/plotting.rst
new file mode 100644
index 000000000..d762a6755
--- /dev/null
+++ b/doc/plotting.rst
@@ -0,0 +1,142 @@
+.. _plotting-module:
+
+*************
+Plotting data
+*************
+
+The Python Control Toolbox contains a number of functions for plotting
+input/output responses in the time and frequency domain, root locus
+diagrams, and other standard charts used in control system analysis.  While
+some legacy functions do both analysis and plotting, the standard pattern
+used in the toolbox is to provide a function that performs the basic
+computation (e.g., time or frequency response) and returns and object
+representing the output data.  A separate plotting function, typically
+ending in `_plot` is then used to plot the data.  The plotting function is
+also available via the `plot()` method of the analysis object, allowing the
+following type of calls::
+
+  step_response(sys).plot()
+  frequency_response(sys).plot()  # implementation pending
+  nyquist_curve(sys).plot()       # implementation pending
+  rootlocus_curve(sys).plot()     # implementation pending
+
+Time response data
+==================
+
+Input/output time responses are produced one of several python-control
+functions: :func:`~control.forced_response`,
+:func:`~control.impulse_response`, :func:`~control.initial_response`,
+:func:`~control.input_output_response`, :func:`~control.step_response`.
+Each of these return a :class:`~control.TimeResponseData` object, which
+contains the time, input, state, and output vectors associated with the
+simulation. Time response data can be plotted with the
+:func:`~control.time_response_plot` function, which is also available as
+the :func:`~control.TimeResponseData.plot` method.  For example, the step
+response for a two-input, two-output can be plotted using the commands::
+
+  sys_mimo = ct.tf2ss(
+      [[[1], [0.1]], [[0.2], [1]]],
+      [[[1, 0.6, 1], [1, 1, 1]], [[1, 0.4, 1], [1, 2, 1]]], name="MIMO")
+  response = step_response(sys)
+  response.plot()
+
+which produces the following plot:
+
+.. image:: timeplot-mimo_step-default.png
+
+The  :class:`~control.TimeResponseData` object can also be used to access
+the data from the simulation::
+
+  time, outputs, inputs = response.time, response.outputs, response.inputs
+  fig, axs = plt.subplots(2, 2)
+  for i in range(2):
+      for j in range(2):
+          axs[i, j].plot(time, outputs[i, j])
+
+A number of options are available in the `plot` method to customize
+the appearance of input output data.  For data produced by the
+:func:`~control.impulse_response` and :func:`~control.step_response`
+commands, the inputs are not shown.  This behavior can be changed
+using the `plot_inputs` keyword.  It is also possible to combine
+multiple lines onto a single graph, using either the `overlay_signals`
+keyword (which puts all outputs out a single graph and all inputs on a
+single graph) or the `overlay_traces` keyword, which puts different
+traces (e.g., corresponding to step inputs in different channels) on
+the same graph, with appropriate labeling via a legend on selected
+axes.
+
+For example, using `plot_input=True` and `overlay_signals=True` yields the
+following plot::
+
+      ct.step_response(sys_mimo).plot(
+        plot_inputs=True, overlay_signals=True,
+        title="Step response for 2x2 MIMO system " + 
+        "[plot_inputs, overlay_signals]")
+
+.. image:: timeplot-mimo_step-pi_cs.png
+
+Input/output response plots created with either the
+:func:`~control.forced_response` or the
+:func:`~control.input_output_response` functions include the input signals by
+default. These can be plotted on separate axes, but also "overlaid" on the
+output axes (useful when the input and output signals are being compared to
+each other).  The following plot shows the use of `plot_inputs='overlay'`
+as well as the ability to reposition the legends using the `legend_map`
+keyword::
+
+    timepts = np.linspace(0, 10, 100)
+    U = np.vstack([np.sin(timepts), np.cos(2*timepts)])
+    ct.input_output_response(sys_mimo, timepts, U).plot(
+        plot_inputs='overlay',
+        legend_map=np.array([['lower right'], ['lower right']]),
+        title="I/O response for 2x2 MIMO system " +
+        "[plot_inputs='overlay', legend_map]")
+  
+.. image:: timeplot-mimo_ioresp-ov_lm.png
+
+Another option that is available is to use the `transpose` keyword so that
+instead of plotting the outputs on the top and inputs on the bottom, the
+inputs are plotted on the left and outputs on the right, as shown in the
+following figure::
+
+    U1 = np.vstack([np.sin(timepts), np.cos(2*timepts)])
+    resp1 = ct.input_output_response(sys_mimo, timepts, U1)
+
+    U2 = np.vstack([np.cos(2*timepts), np.sin(timepts)])
+    resp2 = ct.input_output_response(sys_mimo, timepts, U2)
+
+    ct.combine_time_responses(
+        [resp1, resp2], trace_labels=["Scenario #1", "Scenario #2"]).plot(
+            transpose=True,
+            title="I/O responses for 2x2 MIMO system, multiple traces "
+            "[transpose]")
+	    
+.. image:: timeplot-mimo_ioresp-mt_tr.png
+
+This figure also illustrates the ability to create "multi-trace" plots
+using the :func:`~control.combine_time_responses` function.  The line
+properties that are used when combining signals and traces are set by
+the `input_props`, `output_props` and `trace_props` parameters for
+:func:`~control.time_response_plot`.
+
+Additional customization is possible using the `input_props`,
+`output_props`, and `trace_props` keywords to set complementary line colors
+and styles for various signals and traces::
+
+    out = ct.step_response(sys_mimo).plot(
+        plot_inputs='overlay', overlay_signals=True, overlay_traces=True,
+        output_props=[{'color': c} for c in ['blue', 'orange']],
+        input_props=[{'color': c} for c in ['red', 'green']],
+        trace_props=[{'linestyle': s} for s in ['-', '--']])
+
+.. image:: timeplot-mimo_step-linestyle.png
+
+Plotting functions
+==================
+
+.. autosummary::
+   :toctree: generated/
+
+   ~control.time_response_plot
+   ~control.combine_time_responses
+   ~control.get_plot_axes
diff --git a/doc/timeplot-mimo_ioresp-mt_tr.png b/doc/timeplot-mimo_ioresp-mt_tr.png
new file mode 100644
index 000000000..e4c800086
Binary files /dev/null and b/doc/timeplot-mimo_ioresp-mt_tr.png differ
diff --git a/doc/timeplot-mimo_ioresp-ov_lm.png b/doc/timeplot-mimo_ioresp-ov_lm.png
new file mode 100644
index 000000000..27dd89159
Binary files /dev/null and b/doc/timeplot-mimo_ioresp-ov_lm.png differ
diff --git a/doc/timeplot-mimo_step-default.png b/doc/timeplot-mimo_step-default.png
new file mode 100644
index 000000000..877764fbf
Binary files /dev/null and b/doc/timeplot-mimo_step-default.png differ
diff --git a/doc/timeplot-mimo_step-linestyle.png b/doc/timeplot-mimo_step-linestyle.png
new file mode 100644
index 000000000..9685ea6fa
Binary files /dev/null and b/doc/timeplot-mimo_step-linestyle.png differ
diff --git a/doc/timeplot-mimo_step-pi_cs.png b/doc/timeplot-mimo_step-pi_cs.png
new file mode 100644
index 000000000..6046c8cce
Binary files /dev/null and b/doc/timeplot-mimo_step-pi_cs.png differ