diff --git a/control/__init__.py b/control/__init__.py
index 340370c8b..cfc23ed19 100644
--- a/control/__init__.py
+++ b/control/__init__.py
@@ -54,10 +54,15 @@
 
 Available subpackages
 ---------------------
-flatsys
-    Differentially flat systems
-optimal
-    Optimization-based control
+
+The main control package includes the most commpon functions used in
+analysis, design, and simulation of feedback control systems.  Several
+additional subpackages are available that provide more specialized
+functionality:
+
+* :mod:`~control.flatsys`: Differentially flat systems
+* :mod:`~control.matlab`: MATLAB compatibility module
+* :mod:`~control.optimal`: Optimization-based control
 
 """
 
diff --git a/control/flatsys/systraj.py b/control/flatsys/systraj.py
index c9bde6d7a..0fbd4e982 100644
--- a/control/flatsys/systraj.py
+++ b/control/flatsys/systraj.py
@@ -40,7 +40,7 @@
 from ..timeresp import TimeResponseData
 
 class SystemTrajectory:
-    """Class representing a system trajectory.
+    """Class representing a trajectory for a flat system.
 
     The `SystemTrajectory` class is used to represent the
     trajectory of a (differentially flat) system.  Used by the
diff --git a/control/iosys.py b/control/iosys.py
index dca00d3e5..a468d79e3 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -1704,6 +1704,15 @@ def input_output_response(
        start with zero initial condition since this can be specified as
        [xsys_0, 0].  A warning is issued if the initial conditions are padded
        and and the final listed initial state is not zero.
+    
+    2. If discontinuous inputs are given, the underlying SciPy numerical
+       integration algorithms can sometimes produce erroneous results due
+       to the default tolerances that are used.  The `ivp_method` and
+       `ivp_keywords` parameters can be used to tune the ODE solver and
+       produce better results.  In particular, using 'LSODA' as the
+       `ivp_method` or setting the `rtol` parameter to a smaller value
+       (e.g. using `ivp_kwargs={'rtol': 1e-4}`) can provide more accurate
+       results.
 
     """
     #
diff --git a/control/optimal.py b/control/optimal.py
index 5dfbb34f0..50145324f 100644
--- a/control/optimal.py
+++ b/control/optimal.py
@@ -3,7 +3,7 @@
 # RMM, 11 Feb 2021
 #
 
-"""The :mod:`~control.optimal` module provides support for optimization-based
+"""The :mod:`control.optimal` module provides support for optimization-based
 controllers for nonlinear systems with state and input constraints.
 
 The docstring examples assume that the following import commands::
diff --git a/control/timeresp.py b/control/timeresp.py
index bd8595cfe..2e25331d1 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -820,7 +820,7 @@ def shape_matches(s_legal, s_actual):
 # Forced response of a linear system
 def forced_response(sys, T=None, U=0., X0=0., transpose=False,
                     interpolate=False, return_x=None, squeeze=None):
-    """Simulate the output of a linear system.
+    """Compute the output of a linear system given the input.
 
     As a convenience for parameters `U`, `X0`:
     Numbers (scalars) are converted to constant arrays with the correct shape.
@@ -1616,7 +1616,7 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None,
 def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
                      transpose=False, return_x=False, squeeze=None):
     # pylint: disable=W0622
-    """Initial condition response of a linear system
+    """Compute the initial condition response for a linear system.
 
     If the system has multiple outputs (MIMO), optionally, one output
     may be selected. If no selection is made for the output, all
diff --git a/doc/classes.rst b/doc/classes.rst
index 87ce457de..8564533b3 100644
--- a/doc/classes.rst
+++ b/doc/classes.rst
@@ -58,3 +58,8 @@ Additional classes
    flatsys.SystemTrajectory
    optimal.OptimalControlProblem
    optimal.OptimalControlResult
+   optimal.OptimalEstimationProblem
+   optimal.OptimalEstimationResult
+
+The use of these classes is described in more detail in the
+:ref:`flatsys-module` module and the :ref:`optimal-module` module
diff --git a/doc/conventions.rst b/doc/conventions.rst
index 5dc2e3d76..7c9c1ec6f 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -11,7 +11,6 @@ way that different types of standard information used by the library.
 Throughout this manual, we assume the `control` package has been
 imported as `ct`.
 
-
 LTI system representation
 =========================
 
@@ -132,11 +131,40 @@ constructor for the desired data type using the original system as the sole
 argument or using the explicit conversion functions :func:`ss2tf` and
 :func:`tf2ss`.
 
+Simulating LTI systems
+======================
+
+A number of functions are available for computing the output (and
+state) response of an LTI systems:
+
+.. autosummary::
+   :toctree: generated/
+   :template: custom-class-template.rst
+
+    initial_response
+    step_response
+    impulse_response
+    forced_response
+
+Each of these functions returns a :class:`TimeResponseData` object
+that contains the data for the time response (described in more detail
+in the next section).
+
+The :func:`forced_response` system is the most general and allows by
+the zero initial state response to be simulated as well as the
+response from a non-zero intial condition.
+
+In addition the :func:`input_output_response` function, which handles
+simulation of nonlinear systems and interconnected systems, can be
+used.  For an LTI system, results are generally more accurate using
+the LTI simulation functions above.  The :func:`input_output_response`
+function is described in more detail in the :ref:`iosys-module` section.
+
 .. currentmodule:: control
 .. _time-series-convention:
 
 Time series data
-================
+----------------
 A variety of functions in the library return time series data: sequences of
 values that change over time.  A common set of conventions is used for
 returning such data: columns represent different points in time, rows are
@@ -262,16 +290,16 @@ Selected variables that can be configured, along with their default values:
 
   * freqplot.dB (False): Bode plot magnitude plotted in dB (otherwise powers
     of 10)
-    
+
   * freqplot.deg (True): Bode plot phase plotted in degrees (otherwise radians)
-    
+
   * freqplot.Hz (False): Bode plot frequency plotted in Hertz (otherwise
     rad/sec)
-    
+
   * freqplot.grid (True): Include grids for magnitude and phase plots
-    
+
   * freqplot.number_of_samples (1000): Number of frequency points in Bode plots
-    
+
   * freqplot.feature_periphery_decade (1.0): How many decades to include in
     the frequency range on both sides of features (poles, zeros).