diff --git a/control/dtime.py b/control/dtime.py
index 9dddd86a3..724eafb76 100644
--- a/control/dtime.py
+++ b/control/dtime.py
@@ -90,11 +90,11 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None,
     sysd : linsys
         Discrete time system, with sampling rate Ts
 
-    Additional Parameters
-    ---------------------
+    Other Parameters
+    ----------------
     inputs : int, list of str or None, optional
         Description of the system inputs.  If not specified, the origional
-        system inputs are used.  See :class:`NamedIOSystem` for more
+        system inputs are used.  See :class:`InputOutputSystem` for more
         information.
     outputs : int, list of str or None, optional
         Description of the system outputs.  Same format as `inputs`.
diff --git a/control/iosys.py b/control/iosys.py
index cbbd4cecc..e3af5d9f7 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -2210,17 +2210,16 @@ def linearize(sys, xeq, ueq=None, t=0, params=None, **kw):
         The linearization of the system, as a :class:`~control.LinearIOSystem`
         object (which is also a :class:`~control.StateSpace` object.
 
-    Additional Parameters
-    ---------------------
+    Other Parameters
+    ----------------
     inputs : int, list of str or None, optional
         Description of the system inputs.  If not specified, the origional
-        system inputs are used.  See :class:`NamedIOSystem` for more
+        system inputs are used.  See :class:`InputOutputSystem` for more
         information.
     outputs : int, list of str or None, optional
         Description of the system outputs.  Same format as `inputs`.
     states : int, list of str, or None, optional
         Description of the system states.  Same format as `inputs`.
-
     """
     if not isinstance(sys, InputOutputSystem):
         raise TypeError("Can only linearize InputOutputSystem types")
@@ -2429,7 +2428,10 @@ def rss(states=1, outputs=1, inputs=1, strictly_proper=False, **kwargs):
 
 
 def drss(*args, **kwargs):
-    """Create a stable, discrete-time, random state space system
+    """
+    drss([states, outputs, inputs, strictly_proper])
+
+    Create a stable, discrete-time, random state space system
 
     Create a stable *discrete time* random state space object.  This
     function calls :func:`rss` using either the `dt` keyword provided by
@@ -2467,7 +2469,7 @@ def ss2io(*args, **kwargs):
 
 # Convert a transfer function into an input/output system (wrapper)
 def tf2io(*args, **kwargs):
-    """tf2io(sys)
+    """tf2io(sys[, ...])
 
     Convert a transfer function into an I/O system
 
@@ -2716,8 +2718,8 @@ def interconnect(syslist, connections=None, inplist=None, outlist=None,
     :func:`~control.summing_block` function and the ability to automatically
     interconnect signals with the same names:
 
-    >>> P = control.tf2io(control.tf(1, [1, 0]), inputs='u', outputs='y')
-    >>> C = control.tf2io(control.tf(10, [1, 1]), inputs='e', outputs='u')
+    >>> P = control.tf(1, [1, 0], inputs='u', outputs='y')
+    >>> C = control.tf(10, [1, 1], inputs='e', outputs='u')
     >>> sumblk = control.summing_junction(inputs=['r', '-y'], output='e')
     >>> T = control.interconnect([P, C, sumblk], inputs='r', outputs='y')
 
diff --git a/control/mateqn.py b/control/mateqn.py
index 23ae1e64e..1cf2e65d9 100644
--- a/control/mateqn.py
+++ b/control/mateqn.py
@@ -88,7 +88,9 @@ def sb03md(n, C, A, U, dico, job='X', fact='N', trana='N', ldwork=None):
 
 
 def lyap(A, Q, C=None, E=None, method=None):
-    """X = lyap(A, Q) solves the continuous-time Lyapunov equation
+    """Solves the continuous-time Lyapunov equation
+
+    X = lyap(A, Q) solves
 
         :math:`A X + X A^T + Q = 0`
 
@@ -217,7 +219,9 @@ def lyap(A, Q, C=None, E=None, method=None):
 
 
 def dlyap(A, Q, C=None, E=None, method=None):
-    """dlyap(A, Q) solves the discrete-time Lyapunov equation
+    """Solves the discrete-time Lyapunov equation
+
+    X = dlyap(A, Q) solves
 
         :math:`A X A^T - X + Q = 0`
 
@@ -348,8 +352,9 @@ def dlyap(A, Q, C=None, E=None, method=None):
 
 def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
          A_s="A", B_s="B", Q_s="Q", R_s="R", S_s="S", E_s="E"):
-    """X, L, G = care(A, B, Q, R=None) solves the continuous-time
-    algebraic Riccati equation
+    """Solves the continuous-time algebraic Riccati equation
+
+    X, L, G = care(A, B, Q, R=None) solves
 
         :math:`A^T X + X A - X B R^{-1} B^T X + Q = 0`
 
@@ -505,9 +510,11 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
 
 def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None,
          A_s="A", B_s="B", Q_s="Q", R_s="R", S_s="S", E_s="E"):
-    """X, L, G = dare(A, B, Q, R) solves the discrete-time algebraic Riccati
+    """Solves the discrete-time algebraic Riccati
     equation
 
+    X, L, G = dare(A, B, Q, R) solves
+
         :math:`A^T X A - X - A^T X B (B^T X B + R)^{-1} B^T X A + Q = 0`
 
     where A and Q are square matrices of the same dimension. Further, Q
diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py
index 80f2a0a65..fc0cd445b 100644
--- a/control/matlab/__init__.py
+++ b/control/matlab/__init__.py
@@ -84,6 +84,7 @@
 from ..rlocus import rlocus
 from ..dtime import c2d
 from ..sisotool import sisotool
+from ..stochsys import lqe, dlqe
 
 # Functions that are renamed in MATLAB
 pole, zero = poles, zeros
diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py
index b1fa24bb0..58b5e589d 100644
--- a/control/matlab/timeresp.py
+++ b/control/matlab/timeresp.py
@@ -7,8 +7,7 @@
 __all__ = ['step', 'stepinfo', 'impulse', 'initial', 'lsim']
 
 def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
-    '''
-    Step response of a linear system
+    '''Step response of a linear system
 
     If the system has multiple inputs or outputs (MIMO), one input has
     to be selected for the simulation.  Optionally, one output may be
@@ -20,19 +19,14 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
     ----------
     sys: StateSpace, or TransferFunction
         LTI system to simulate
-
     T: array-like or number, optional
         Time vector, or simulation time duration if a number (time vector is
         autocomputed if not given)
-
     X0: array-like or number, optional
         Initial condition (default = 0)
-
         Numbers are converted to constant arrays with the correct shape.
-
     input: int
         Index of the input that will be used in this simulation.
-
     output: int
         If given, index of the output that is returned by this simulation.
 
@@ -40,15 +34,12 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
     -------
     yout: array
         Response of the system
-
     T: array
         Time values of the output
-
     xout: array (if selected)
         Individual response of each x variable
 
 
-
     See Also
     --------
     lsim, initial, impulse
@@ -67,8 +58,7 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
 
 def stepinfo(sysdata, T=None, yfinal=None, SettlingTimeThreshold=0.02,
              RiseTimeLimits=(0.1, 0.9)):
-    """
-    Step response characteristics (Rise time, Settling Time, Peak and others).
+    """Step response characteristics (Rise time, Settling Time, Peak and others)
 
     Parameters
     ----------
@@ -137,8 +127,7 @@ def stepinfo(sysdata, T=None, yfinal=None, SettlingTimeThreshold=0.02,
     return S
 
 def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
-    '''
-    Impulse response of a linear system
+    '''Impulse response of a linear system
 
     If the system has multiple inputs or outputs (MIMO), one input has
     to be selected for the simulation.  Optionally, one output may be
@@ -150,19 +139,15 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
     ----------
     sys: StateSpace, TransferFunction
         LTI system to simulate
-
     T: array-like or number, optional
         Time vector, or simulation time duration if a number (time vector is
         autocomputed if not given)
-
     X0: array-like or number, optional
         Initial condition (default = 0)
 
         Numbers are converted to constant arrays with the correct shape.
-
     input: int
         Index of the input that will be used in this simulation.
-
     output: int
         Index of the output that will be used in this simulation.
 
@@ -170,10 +155,8 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
     -------
     yout: array
         Response of the system
-
     T: array
         Time values of the output
-
     xout: array (if selected)
         Individual response of each x variable
 
@@ -193,8 +176,7 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
     return (out[1], out[0], out[2]) if return_x else (out[1], out[0])
 
 def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
-    '''
-    Initial condition response of a linear system
+    '''Initial condition response of a linear system
 
     If the system has multiple outputs (?IMO), optionally, one output
     may be selected. If no selection is made for the output, all
@@ -204,20 +186,16 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
     ----------
     sys: StateSpace, or TransferFunction
         LTI system to simulate
-
     T: array-like or number, optional
         Time vector, or simulation time duration if a number (time vector is
         autocomputed if not given)
-
     X0: array-like object or number, optional
         Initial condition (default = 0)
 
         Numbers are converted to constant arrays with the correct shape.
-
     input: int
         This input is ignored, but present for compatibility with step
         and impulse.
-
     output: int
         If given, index of the output that is returned by this simulation.
 
@@ -225,10 +203,8 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
     -------
     yout: array
         Response of the system
-
     T: array
         Time values of the output
-
     xout: array (if selected)
         Individual response of each x variable
 
@@ -250,8 +226,7 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
 
 
 def lsim(sys, U=0., T=None, X0=0.):
-    '''
-    Simulate the output of a linear system.
+    '''Simulate the output of a linear system
 
     As a convenience for parameters `U`, `X0`:
     Numbers (scalars) are converted to constant arrays with the correct shape.
@@ -261,16 +236,13 @@ def lsim(sys, U=0., T=None, X0=0.):
     ----------
     sys: LTI (StateSpace, or TransferFunction)
         LTI system to simulate
-
     U: array-like or number, optional
         Input array giving input at each time `T` (default = 0).
 
         If `U` is ``None`` or ``0``, a special algorithm is used. This special
         algorithm is faster than the general algorithm, which is used otherwise.
-
     T: array-like, optional for discrete LTI `sys`
         Time steps at which the input is defined; values must be evenly spaced.
-
     X0: array-like or number, optional
         Initial condition (default = 0).
 
diff --git a/control/matlab/wrappers.py b/control/matlab/wrappers.py
index 8eafdaad2..6a6c4ad35 100644
--- a/control/matlab/wrappers.py
+++ b/control/matlab/wrappers.py
@@ -178,8 +178,7 @@ def ngrid():
 
 
 def dcgain(*args):
-    '''
-    Compute the gain of the system in steady state.
+    '''Compute the gain of the system in steady state
 
     The function takes either 1, 2, 3, or 4 parameters:
 
diff --git a/control/sisotool.py b/control/sisotool.py
index 781fabf40..598bdcf09 100644
--- a/control/sisotool.py
+++ b/control/sisotool.py
@@ -209,10 +209,6 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
 
     C_f = Kp + Ki/s + Kd*s/(tau*s + 1).
 
-    If `plant` is a discrete-time system, then the proportional, integral, and
-    derivative terms are given instead by Kp, Ki*dt/2*(z+1)/(z-1), and
-    Kd/dt*(z-1)/z, respectively.
-
     ::
 
           ------> C_ff ------    d
@@ -224,6 +220,10 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
               |             ----- C_b <-------|
               ---------------------------------
 
+    If `plant` is a discrete-time system, then the proportional, integral, and
+    derivative terms are given instead by Kp, Ki*dt/2*(z+1)/(z-1), and
+    Kd/dt*(z-1)/z, respectively.
+
     It is also possible to move the derivative term into the feedback path
     `C_b` using `derivative_in_feedback_path=True`. This may be desired to
     avoid that the plant is subject to an impulse function when the reference
diff --git a/control/statesp.py b/control/statesp.py
index 7843cb33f..0303ff70a 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -807,7 +807,7 @@ def __rdiv__(self, other):
             "StateSpace.__rdiv__ is not implemented yet.")
 
     def __call__(self, x, squeeze=None, warn_infinite=True):
-        """Evaluate system's transfer function at complex frequency.
+        """Evaluate system's frequency response at complex frequencies.
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
         continuous-time systems and `z` for discrete-time systems.
@@ -1344,8 +1344,8 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None,
         sysd : StateSpace
             Discrete-time system, with sampling rate Ts
 
-        Additional Parameters
-        ---------------------
+        Other Parameters
+        ----------------
         inputs : int, list of str or None, optional
             Description of the system inputs.  If not specified, the origional
             system inputs are used.  See :class:`InputOutputSystem` for more
@@ -1520,6 +1520,7 @@ def output(self, t, x, u=None, params=None):
 
 
 # TODO: add discrete time check
+# TODO: copy signal names
 def _convert_to_statespace(sys):
     """Convert a system to state space form (if needed).
 
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 84188a63f..490132952 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1136,11 +1136,11 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None,
         sysd : TransferFunction system
             Discrete-time system, with sample period Ts
 
-        Additional Parameters
-        ---------------------
+        Other Parameters
+        ----------------
         inputs : int, list of str or None, optional
             Description of the system inputs.  If not specified, the origional
-            system inputs are used.  See :class:`NamedIOSystem` for more
+            system inputs are used.  See :class:`InputOutputSystem` for more
             information.
         outputs : int, list of str or None, optional
             Description of the system outputs.  Same format as `inputs`.
@@ -1585,7 +1585,7 @@ def tf(*args, **kwargs):
     else:
         raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
-
+# TODO: copy signal names
 def ss2tf(*args, **kwargs):
     """ss2tf(sys)
 
diff --git a/doc/control.rst b/doc/control.rst
index 172790f83..54e233746 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -20,6 +20,8 @@ System creation
     frd
     rss
     drss
+    NonlinearIOSystem
+
 
 System interconnections
 =======================
@@ -29,12 +31,11 @@ System interconnections
     append
     connect
     feedback
+    interconnect
     negate
     parallel
     series
 
-See also the :ref:`iosys-module` module, which can be used to create and
-interconnect nonlinear input/output systems.
 
 Frequency domain plotting
 =========================
@@ -78,8 +79,7 @@ Control system analysis
 
     dcgain
     describing_function
-    evalfr
-    freqresp
+    frequency_response
     get_input_ff_index
     get_output_fb_index
     ispassive
@@ -91,6 +91,8 @@ Control system analysis
     pzmap
     root_locus
     sisotool
+    StateSpace.__call__
+    TransferFunction.__call__
 
 
 
@@ -100,10 +102,10 @@ Matrix computations
    :toctree: generated/
 
     care
+    ctrb
     dare
-    lyap
     dlyap
-    ctrb
+    lyap
     obsv
     gram
 
@@ -141,7 +143,6 @@ Nonlinear system support
 
     describing_function
     find_eqpt
-    interconnect
     linearize
     input_output_response
     ss2io
@@ -193,3 +194,5 @@ Utility functions and conversions
     use_fbs_defaults
     use_matlab_defaults
     use_numpy_matrix
+
+
diff --git a/doc/matlab.rst b/doc/matlab.rst
index c14a67e1f..eac1d157a 100644
--- a/doc/matlab.rst
+++ b/doc/matlab.rst
@@ -86,6 +86,9 @@ Compensator design
    sisotool
    place
    lqr
+   dlqr
+   lqe
+   dlqe
 
 State-space (SS) models
 =======================