Repository: python-control/python-control
Pull request #920 (2023-07-03 01:13:16)

Time response plots
This PR implements time response plots for the various simulation routines in python-control via a new function `time_response_plot()` and via the `plot()` method on the `TimeResponseData` class.  This is aligned with the new plotting paradigm in #645.

**Examples** (from the user documentation):

Input/output time responses are produced using one of several python-control functions: `forced_response`, `impulse_response`, `initial_response()`, `input_output_response()`, and `step_response()`. Each of these return a `TimeResponseData` object, which contains the time, input, state, and output vectors associated with the simulation. Time response data can be plotted with the `time_response_plot()` function, which is also available as the `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:

![timeplot-mimo_step-default](https://github.com/python-control/python-control/assets/293362/18070852-c873-418d-a9bd-61362b7cfacf)

The `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 `impulse_response()` and `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]")
```
![timeplot-mimo_step-pi_cs](https://github.com/python-control/python-control/assets/293362/8ea7216f-209a-4ced-98e8-b246abf48ed6)

Input/output response plots created with either the `forced_response()` or the `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]")
```
![timeplot-mimo_ioresp-ov_lm](https://github.com/python-control/python-control/assets/293362/04e35e16-4b24-4cd3-8aca-72d4618e2756)

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_traces(
    [resp1, resp2], trace_labels=["Scenario #1", "Scenario #2"]).plot(
        transpose=True,
        title="I/O responses for 2x2 MIMO system, multiple traces "
        "[transpose]")
```
![timeplot-mimo_ioresp-mt_tr](https://github.com/python-control/python-control/assets/293362/0c3941e6-d88b-46f7-b779-305df546ebe3)

This figure also illustrates the ability to create “multi-trace” plots using the `combine_traces()` function.

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 ['-', '--']])
```
![timeplot-mimo_step-linestyle](https://github.com/python-control/python-control/assets/293362/f054d72a-d629-4aa4-96da-9632f24b2c5c)


coveralls: 
[![Coverage Status](https://coveralls.io/builds/61092978/badge)](https://coveralls.io/builds/61092978)

coverage: 94.967% (+0.1%) from 94.849% when pulling **dd6983ef83e5cf20e6a9e7ed278d377e5594f1ff on murrayrm:time_plots-20Jun2023** into **d1530357f1d57837756ea7adcb0ac592d473eb93 on python-control:main**.


Patch:
diff --git a/control/__init__.py b/control/__init__.py
index 9781ab80..120d1632 100644
--- a/control/__init__.py
+++ b/control/__init__.py
@@ -77,6 +77,10 @@ from .statesp import *
 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 .pzmap import *
 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 987693c2..1ed8b5dd 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 c540decb..f11d08a3 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -1488,6 +1488,7 @@ def input_output_response(
         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 input_output_response(
     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/timeplot.py b/control/timeplot.py
new file mode 100644
index 00000000..6409a666
--- /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 defdbdf4..bc13ad0d 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -81,6 +81,7 @@ from copy import copy
 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 @@ class TimeResponseData:
             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 @@ class TimeResponseData:
             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 @@ class TimeResponseData:
             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 @@ class TimeResponseData:
         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 @@ class TimeResponseData:
         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 @@ class TimeResponseData:
         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 @@ class TimeResponseData:
 
         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 6e101234..88a1b7ba 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 98b18428..ec556e7c 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 00000000..d762a675
--- /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


Test patch:
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index 1abc7547..19f4bb62 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -26,6 +26,7 @@ import control.tests.optimal_test as optimal_test
 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 @@ kwarg_unittest = {
     '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 @@ kwarg_unittest = {
     '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 00000000..7cdde5c5
--- /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')