Repository: python-control/python-control

Issue:
step_info error for MIMO system 
I try with this MIMO system, the function broke at line 614. The step_info function convert the system to SIMO, but cannot calculate the step response characteristic parameters.

```
from control import ss step_response step_info

A = [[0.68, -0.34], [0.34, 0.68]];
B = [[0.18, -0.05], [0.04, 0.11]];
C = [[0, -1.53], [-1.12, -1.10]];
D = [[0, 0], [0.06, -0.37]];
m_sys = ss(A,B,C,D,0.2);
mt,my=step_response(m_sys)

step_info(m_sys)

```


>   File ".../timeresp.py", line 614, in step_info
>     if((yout[i] <= inf_margin) | (yout[i] >= sup_margin)):
> 
>  IndexError: index 24 is out of bounds for axis 0 with size 2
> 

Discussion:

Patch:
diff --git a/control/iosys.py b/control/iosys.py
index e28de59f..7ed4c8b0 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -358,88 +358,26 @@ class InputOutputSystem(object):
         if (warning):
             warn("Parameters passed to InputOutputSystem ignored.")
 
-    def _rhs(self, t, x, u, params={}):
+    def _rhs(self, t, x, u):
         """Evaluate right hand side of a differential or difference equation.
 
         Private function used to compute the right hand side of an
-        input/output system model. Intended for fast
-        evaluation; for a more user-friendly interface
-        you may want to use :meth:`dynamics`.
+        input/output system model.
 
         """
         NotImplemented("Evaluation not implemented for system of type ",
                        type(self))
 
-    def dynamics(self, t, x, u):
-        """Compute the dynamics of a differential or difference equation.
-
-        Given time `t`, input `u` and state `x`, returns the value of the
-        right hand side of the dynamical system. If the system is continuous,
-        returns the time derivative
-
-            dx/dt = f(t, x, u)
-
-        where `f` is the system's (possibly nonlinear) dynamics function.
-        If the system is discrete-time, returns the next value of `x`:
-
-            x[t+dt] = f(t, x[t], u[t])
-
-        Where `t` is a scalar.
-
-        The inputs `x` and `u` must be of the correct length.
-
-        Parameters
-        ----------
-        t : float
-            the time at which to evaluate
-        x : array_like
-            current state
-        u : array_like
-            input
-
-        Returns
-        -------
-        dx/dt or x[t+dt] : ndarray
-        """
-        return self._rhs(t, x, u)
-
     def _out(self, t, x, u, params={}):
         """Evaluate the output of a system at a given state, input, and time
 
         Private function used to compute the output of of an input/output
-        system model given the state, input, parameters. Intended for fast
-        evaluation; for a more user-friendly interface you may want to use
-        :meth:`output`.
+        system model given the state, input, parameters, and time.
 
         """
         # If no output function was defined in subclass, return state
         return x
 
-    def output(self, t, x, u):
-        """Compute the output of the system
-
-        Given time `t`, input `u` and state `x`, returns the output of the
-        system:
-
-            y = g(t, x, u)
-
-        The inputs `x` and `u` must be of the correct length.
-
-        Parameters
-        ----------
-        t : float
-            the time at which to evaluate
-        x : array_like
-            current state
-        u : array_like
-            input
-
-        Returns
-        -------
-        y : ndarray
-        """
-        return self._out(t, x, u)
-
     def set_inputs(self, inputs, prefix='u'):
         """Set the number/names of the system inputs.
 
diff --git a/control/margins.py b/control/margins.py
index 0b53f26e..d22b6a0f 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -218,7 +218,7 @@ def _likely_numerical_inaccuracy(sys):
         # * z**(-p_q)
         x = [1] + [0] * (-p_q)
         p1 = np.polymul(p1, x)
-    return np.linalg.norm(p1) < 1e-4 * np.linalg.norm(p2)
+    return np.linalg.norm(p1) < 1e-3 * np.linalg.norm(p2)
 
 # Took the framework for the old function by
 # Sawyer B. Fuller <minster@uw.edu>, removed a lot of the innards
@@ -337,8 +337,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0, method='best'):
         if isinstance(sys, xferfcn.TransferFunction) and not sys.isctime():
             if _likely_numerical_inaccuracy(sys):
                 warn("stability_margins: Falling back to 'frd' method "
-                "because of chance of numerical inaccuracy in 'poly' method.",
-                stacklevel=2)
+                "because of chance of numerical inaccuracy in 'poly' method.")
                 omega_sys = freqplot._default_frequency_range(sys)
                 omega_sys = omega_sys[omega_sys < np.pi / sys.dt]
                 sys = frdata.FRD(sys, omega_sys, smooth=True)
diff --git a/control/mateqn.py b/control/mateqn.py
index 28b01d28..0b129fd9 100644
--- a/control/mateqn.py
+++ b/control/mateqn.py
@@ -5,6 +5,9 @@
 #
 # Author: Bjorn Olofsson
 
+# Python 3 compatibility (needs to go here)
+from __future__ import print_function
+
 # Copyright (c) 2011, All rights reserved.
 
 # Redistribution and use in source and binary forms, with or without
@@ -41,34 +44,6 @@ from scipy.linalg import eigvals, solve_discrete_are, solve
 from .exception import ControlSlycot, ControlArgument
 from .statesp import _ssmatrix
 
-# Make sure we have access to the right slycot routines
-try:
-    from slycot import sb03md57
-    # wrap without the deprecation warning
-    def sb03md(n, C, A, U, dico, job='X',fact='N',trana='N',ldwork=None):
-        ret = sb03md57(A, U, C, dico, job, fact, trana, ldwork)
-        return ret[2:]
-except ImportError:
-    try:
-        from slycot import sb03md
-    except ImportError:
-        sb03md = None
-
-try:
-    from slycot import sb04md
-except ImportError:
-    sb04md = None
-
-try:
-    from slycot import sb04qd
-except ImportError:
-    sb0qmd = None
-
-try:
-    from slycot import sg03ad
-except ImportError:
-    sb04ad = None
-
 __all__ = ['lyap', 'dlyap', 'dare', 'care']
 
 #
@@ -118,11 +93,16 @@ def lyap(A, Q, C=None, E=None):
     state space operations.  See :func:`~control.use_numpy_matrix`.
     """
 
-    if sb03md is None:
+    # Make sure we have access to the right slycot routines
+    try:
+        from slycot import sb03md
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sb03md'")
-    if sb04md is None:
-        raise ControlSlycot("can't find slycot module 'sb04md'")
 
+    try:
+        from slycot import sb04md
+    except ImportError:
+        raise ControlSlycot("can't find slycot module 'sb04md'")
 
     # Reshape 1-d arrays
     if len(shape(A)) == 1:
@@ -299,11 +279,19 @@ def dlyap(A, Q, C=None, E=None):
     of the same dimension. """
 
     # Make sure we have access to the right slycot routines
-    if sb03md is None:
+    try:
+        from slycot import sb03md
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sb03md'")
-    if sb04qd is None:
+
+    try:
+        from slycot import sb04qd
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sb04qd'")
-    if sg03ad is None:
+
+    try:
+        from slycot import sg03ad
+    except ImportError:
         raise ControlSlycot("can't find slycot module 'sg03ad'")
 
     # Reshape 1-d arrays
diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py
index 31b761bc..b1fa24bb 100644
--- a/control/matlab/timeresp.py
+++ b/control/matlab/timeresp.py
@@ -64,38 +64,59 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
                         transpose=True, return_x=return_x)
     return (out[1], out[0], out[2]) if return_x else (out[1], out[0])
 
-def stepinfo(sys, T=None, SettlingTimeThreshold=0.02,
+
+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).
 
     Parameters
     ----------
-    sys: StateSpace, or TransferFunction
-        LTI system to simulate
-
-    T: array-like or number, optional
+    sysdata : StateSpace or TransferFunction or array_like
+        The system data. Either LTI system to similate (StateSpace,
+        TransferFunction), or a time series of step response data.
+    T : array_like or float, optional
         Time vector, or simulation time duration if a number (time vector is
-        autocomputed if not given)
-
-    SettlingTimeThreshold: float value, optional
+        autocomputed if not given).
+        Required, if sysdata is a time series of response data.
+    yfinal : scalar or array_like, optional
+        Steady-state response. If not given, sysdata.dcgain() is used for
+        systems to simulate and the last value of the the response data is
+        used for a given time series of response data. Scalar for SISO,
+        (noutputs, ninputs) array_like for MIMO systems.
+    SettlingTimeThreshold : float, optional
         Defines the error to compute settling time (default = 0.02)
-
-    RiseTimeLimits: tuple (lower_threshold, upper_theshold)
+    RiseTimeLimits : tuple (lower_threshold, upper_theshold)
         Defines the lower and upper threshold for RiseTime computation
 
     Returns
     -------
-    S: a dictionary containing:
-        RiseTime: Time from 10% to 90% of the steady-state value.
-        SettlingTime: Time to enter inside a default error of 2%
-        SettlingMin: Minimum value after RiseTime
-        SettlingMax: Maximum value after RiseTime
-        Overshoot: Percentage of the Peak relative to steady value
-        Undershoot: Percentage of undershoot
-        Peak: Absolute peak value
-        PeakTime: time of the Peak
-        SteadyStateValue: Steady-state value
+    S : dict or list of list of dict
+        If `sysdata` corresponds to a SISO system, S is a dictionary
+        containing:
+
+        RiseTime:
+            Time from 10% to 90% of the steady-state value.
+        SettlingTime:
+            Time to enter inside a default error of 2%
+        SettlingMin:
+            Minimum value after RiseTime
+        SettlingMax:
+            Maximum value after RiseTime
+        Overshoot:
+            Percentage of the Peak relative to steady value
+        Undershoot:
+            Percentage of undershoot
+        Peak:
+            Absolute peak value
+        PeakTime:
+            time of the Peak
+        SteadyStateValue:
+            Steady-state value
+
+        If `sysdata` corresponds to a MIMO system, `S` is a 2D list of dicts.
+        To get the step response characteristics from the j-th input to the
+        i-th output, access ``S[i][j]``
 
 
     See Also
@@ -105,11 +126,13 @@ def stepinfo(sys, T=None, SettlingTimeThreshold=0.02,
     Examples
     --------
     >>> S = stepinfo(sys, T)
-    '''
+    """
     from ..timeresp import step_info
 
     # Call step_info with MATLAB defaults
-    S = step_info(sys, T, None, SettlingTimeThreshold, RiseTimeLimits)
+    S = step_info(sysdata, T=T, T_num=None, yfinal=yfinal,
+                  SettlingTimeThreshold=SettlingTimeThreshold,
+                  RiseTimeLimits=RiseTimeLimits)
 
     return S
 
diff --git a/control/statefbk.py b/control/statefbk.py
index 0017412a..7106449a 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -47,25 +47,6 @@ from .mateqn import care
 from .statesp import _ssmatrix
 from .exception import ControlSlycot, ControlArgument, ControlDimension
 
-# Make sure we have access to the right slycot routines
-try:
-    from slycot import sb03md57
-    # wrap without the deprecation warning
-    def sb03md(n, C, A, U, dico, job='X',fact='N',trana='N',ldwork=None):
-        ret = sb03md57(A, U, C, dico, job, fact, trana, ldwork)
-        return ret[2:]
-except ImportError:
-    try:
-        from slycot import sb03md
-    except ImportError:
-        sb03md = None
-
-try:
-    from slycot import sb03od
-except ImportError:
-    sb03od = None
-
-
 __all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'lqe',
            'acker']
 
@@ -593,7 +574,7 @@ def gram(sys, type):
         * if `type` is not 'c', 'o', 'cf' or 'of'
         * if system is unstable (sys.A has eigenvalues not in left half plane)
 
-    ControlSlycot
+    ImportError
         if slycot routine sb03md cannot be found
         if slycot routine sb03od cannot be found
 
@@ -633,7 +614,9 @@ def gram(sys, type):
     if type == 'c' or type == 'o':
         # Compute Gramian by the Slycot routine sb03md
         # make sure Slycot is installed
-        if sb03md is None:
+        try:
+            from slycot import sb03md
+        except ImportError:
             raise ControlSlycot("can't find slycot module 'sb03md'")
         if type == 'c':
             tra = 'T'
@@ -651,7 +634,9 @@ def gram(sys, type):
 
     elif type == 'cf' or type == 'of':
         # Compute cholesky factored gramian from slycot routine sb03od
-        if sb03od is None:
+        try:
+            from slycot import sb03od
+        except ImportError:
             raise ControlSlycot("can't find slycot module 'sb03od'")
         tra = 'N'
         n = sys.nstates
diff --git a/control/statesp.py b/control/statesp.py
index 758b91ed..d2b61302 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -1227,101 +1227,12 @@ class StateSpace(LTI):
         return self(0, warn_infinite=warn_infinite) if self.isctime() \
             else self(1, warn_infinite=warn_infinite)
 
-    def dynamics(self, t, x, u=0):
-        """Compute the dynamics of the system
-
-        Given input `u` and state `x`, returns the dynamics of the state-space
-        system. If the system is continuous, returns the time derivative dx/dt
-
-            dx/dt = A x + B u
-
-        where A and B are the state-space matrices of the system. If the
-        system is discrete-time, returns the next value of `x`:
-
-            x[t+dt] = A x[t] + B u[t]
-
-        The inputs `x` and `u` must be of the correct length for the system.
-
-        The first argument `t` is ignored because :class:`StateSpace` systems
-        are time-invariant. It is included so that the dynamics can be passed
-        to most numerical integrators, such as :func:`scipy.integrate.solve_ivp`
-        and for consistency with :class:`IOSystem` systems.
-
-        Parameters
-        ----------
-        t : float (ignored)
-            time
-        x : array_like
-            current state
-        u : array_like (optional)
-            input, zero if omitted
-
-        Returns
-        -------
-        dx/dt or x[t+dt] : ndarray
-        """
-        x = np.reshape(x, (-1, 1)) # force to a column in case matrix
-        if np.size(x) != self.nstates:
-            raise ValueError("len(x) must be equal to number of states")
-        if u is 0:
-            return self.A.dot(x).reshape((-1,)) # return as row vector
-        else: # received t, x, and u, ignore t
-            u = np.reshape(u, (-1, 1)) # force to a column in case matrix
-            if np.size(u) != self.ninputs:
-                raise ValueError("len(u) must be equal to number of inputs")
-            return self.A.dot(x).reshape((-1,)) \
-                 + self.B.dot(u).reshape((-1,)) # return as row vector
-
-    def output(self, t, x, u=0):
-        """Compute the output of the system
-
-        Given input `u` and state `x`, returns the output `y` of the
-        state-space system:
-
-            y = C x + D u
-
-        where A and B are the state-space matrices of the system.
-
-        The first argument `t` is ignored because :class:`StateSpace` systems
-        are time-invariant. It is included so that the dynamics can be passed
-        to most numerical integrators, such as scipy's `integrate.solve_ivp` and
-        for consistency with :class:`IOSystem` systems.
-
-        The inputs `x` and `u` must be of the correct length for the system.
-
-        Parameters
-        ----------
-        t : float (ignored)
-            time
-        x : array_like
-            current state
-        u : array_like (optional)
-            input (zero if omitted)
-
-        Returns
-        -------
-        y : ndarray
-        """
-        x = np.reshape(x, (-1, 1)) # force to a column in case matrix
-        if np.size(x) != self.nstates:
-            raise ValueError("len(x) must be equal to number of states")
-
-        if u is 0:
-            return self.C.dot(x).reshape((-1,)) # return as row vector
-        else: # received t, x, and u, ignore t
-            u = np.reshape(u, (-1, 1)) # force to a column in case matrix
-            if np.size(u) != self.ninputs:
-                raise ValueError("len(u) must be equal to number of inputs")
-            return self.C.dot(x).reshape((-1,)) \
-                 + self.D.dot(u).reshape((-1,)) # return as row vector
-
     def _isstatic(self):
         """True if and only if the system has no dynamics, that is,
         if A and B are zero. """
         return not np.any(self.A) and not np.any(self.B)
 
 
-
 # TODO: add discrete time check
 def _convert_to_statespace(sys, **kw):
     """Convert a system to state space form (if needed).
diff --git a/control/timeresp.py b/control/timeresp.py
index 9c7b7a99..630eff03 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -1,9 +1,7 @@
-# timeresp.py - time-domain simulation routines
-#
-# This file contains a collection of functions that calculate time
-# responses for linear systems.
+"""
+timeresp.py - time-domain simulation routines.
 
-"""The :mod:`~control.timeresp` module contains a collection of
+The :mod:`~control.timeresp` module contains a collection of
 functions that are used to compute time-domain simulations of LTI
 systems.
 
@@ -21,9 +19,7 @@ evaluated, `U` is a vector of inputs (one for each time point) and
 See :ref:`time-series-convention` for more information on how time
 series data are represented.
 
-"""
-
-"""Copyright (c) 2011 by California Institute of Technology
+Copyright (c) 2011 by California Institute of Technology
 All rights reserved.
 
 Copyright (c) 2011 by Eike Welk
@@ -71,18 +67,17 @@ Date: August 17, 2020
 $Id$
 """
 
-# Libraries that we make use of
-import scipy as sp              # SciPy library (used all over)
-import numpy as np              # NumPy library
-from scipy.linalg import eig, eigvals, matrix_balance, norm
-from numpy import (einsum, maximum, minimum,
-                   atleast_1d)
 import warnings
-from .lti import LTI     # base class of StateSpace, TransferFunction
-from .xferfcn import TransferFunction
-from .statesp import _convert_to_statespace, _mimo2simo, _mimo2siso, ssdata
-from .lti import isdtime, isctime
+
+import numpy as np
+import scipy as sp
+from numpy import einsum, maximum, minimum
+from scipy.linalg import eig, eigvals, matrix_balance, norm
+
 from . import config
+from .lti import isctime, isdtime
+from .statesp import StateSpace, _convert_to_statespace, _mimo2simo, _mimo2siso
+from .xferfcn import TransferFunction
 
 __all__ = ['forced_response', 'step_response', 'step_info', 'initial_response',
            'impulse_response']
@@ -210,7 +205,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     Parameters
     ----------
-    sys : LTI (StateSpace or TransferFunction)
+    sys : StateSpace or TransferFunction
         LTI system to simulate
 
     T : array_like, optional for discrete LTI `sys`
@@ -285,9 +280,9 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     See :ref:`time-series-convention`.
 
     """
-    if not isinstance(sys, LTI):
-        raise TypeError('Parameter ``sys``: must be a ``LTI`` object. '
-                        '(For example ``StateSpace`` or ``TransferFunction``)')
+    if not isinstance(sys, (StateSpace, TransferFunction)):
+        raise TypeError('Parameter ``sys``: must be a ``StateSpace`` or'
+                        ' ``TransferFunction``)')
 
     # If return_x was not specified, figure out the default
     if return_x is None:
@@ -739,43 +734,62 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         squeeze=squeeze, input=input, output=output)
 
 
-def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
-              RiseTimeLimits=(0.1, 0.9)):
-    '''
+def step_info(sysdata, T=None, T_num=None, yfinal=None,
+              SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)):
+    """
     Step response characteristics (Rise time, Settling Time, Peak and others).
 
     Parameters
     ----------
-    sys : SISO dynamic system model. Dynamic systems that you can use include:
-        StateSpace or TransferFunction
-        LTI system to simulate
-
+    sysdata : StateSpace or TransferFunction or array_like
+        The system data. Either LTI system to similate (StateSpace,
+        TransferFunction), or a time series of step response data.
     T : array_like or float, optional
         Time vector, or simulation time duration if a number (time vector is
-        autocomputed if not given, see :func:`step_response` for more detail)
-
+        autocomputed if not given, see :func:`step_response` for more detail).
+        Required, if sysdata is a time series of response data.
     T_num : int, optional
         Number of time steps to use in simulation if T is not provided as an
-        array (autocomputed if not given); ignored if sys is discrete-time.
-
-    SettlingTimeThreshold : float value, optional
+        array; autocomputed if not given; ignored if sysdata is a
+        discrete-time system or a time series or response data.
+    yfinal : scalar or array_like, optional
+        Steady-state response. If not given, sysdata.dcgain() is used for
+        systems to simulate and the last value of the the response data is
+        used for a given time series of response data. Scalar for SISO,
+        (noutputs, ninputs) array_like for MIMO systems.
+    SettlingTimeThreshold : float, optional
         Defines the error to compute settling time (default = 0.02)
-
     RiseTimeLimits : tuple (lower_threshold, upper_theshold)
         Defines the lower and upper threshold for RiseTime computation
 
     Returns
     -------
-    S: a dictionary containing:
-        RiseTime: Time from 10% to 90% of the steady-state value.
-        SettlingTime: Time to enter inside a default error of 2%
-        SettlingMin: Minimum value after RiseTime
-        SettlingMax: Maximum value after RiseTime
-        Overshoot: Percentage of the Peak relative to steady value
-        Undershoot: Percentage of undershoot
-        Peak: Absolute peak value
-        PeakTime: time of the Peak
-        SteadyStateValue: Steady-state value
+    S : dict or list of list of dict
+        If `sysdata` corresponds to a SISO system, S is a dictionary
+        containing:
+
+        RiseTime:
+            Time from 10% to 90% of the steady-state value.
+        SettlingTime:
+            Time to enter inside a default error of 2%
+        SettlingMin:
+            Minimum value after RiseTime
+        SettlingMax:
+            Maximum value after RiseTime
+        Overshoot:
+            Percentage of the Peak relative to steady value
+        Undershoot:
+            Percentage of undershoot
+        Peak:
+            Absolute peak value
+        PeakTime:
+            time of the Peak
+        SteadyStateValue:
+            Steady-state value
+
+        If `sysdata` corresponds to a MIMO system, `S` is a 2D list of dicts.
+        To get the step response characteristics from the j-th input to the
+        i-th output, access ``S[i][j]``
 
 
     See Also
@@ -784,83 +798,163 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
 
     Examples
     --------
-    >>> info = step_info(sys, T)
-    '''
-    
-    if not sys.issiso():
-        sys = _mimo2siso(sys,0,0)
-        warnings.warn(" Internal conversion from a MIMO system to a SISO system,"
-                      " the first input and the first output were used (u1 -> y1);"
-                      " it may not be the result you are looking for")
-
-    if T is None or np.asarray(T).size == 1:
-        T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
-
-    T, yout = step_response(sys, T)
-
-    # Steady state value
-    InfValue = sys.dcgain().real
-
-     # TODO: this could be a function step_info4data(t,y,yfinal)
-    rise_time: float = np.NaN
-    settling_time: float = np.NaN
-    settling_min: float = np.NaN
-    settling_max: float = np.NaN
-    peak_value: float = np.Inf
-    peak_time: float = np.Inf
-    undershoot: float = np.NaN
-    overshoot: float = np.NaN
-    steady_state_value: float = np.NaN
- 
-    if not np.isnan(InfValue) and not np.isinf(InfValue):
-        # SteadyStateValue
-        steady_state_value = InfValue
-        # Peak
-        peak_index = np.abs(yout).argmax()
-        peak_value = np.abs(yout[peak_index])
-        peak_time = T[peak_index]
-
-        sup_margin = (1. + SettlingTimeThreshold) * InfValue
-        inf_margin = (1. - SettlingTimeThreshold) * InfValue
-
-        # RiseTime
-        tr_lower_index = (np.where(np.sign(InfValue.real) * (yout- RiseTimeLimits[0] * InfValue) >= 0 )[0])[0]
-        tr_upper_index = (np.where(np.sign(InfValue.real) * yout >= np.sign(InfValue.real) * RiseTimeLimits[1] * InfValue)[0])[0]
-
-        # SettlingTime
-        settling_time = T[np.where(np.abs(yout-InfValue) >= np.abs(SettlingTimeThreshold*InfValue))[0][-1]+1]
-        # Overshoot and Undershoot
-        y_os = (np.sign(InfValue.real)*yout).max()
-        dy_os = np.abs(y_os) - np.abs(InfValue)
-        if dy_os > 0:
-            overshoot = np.abs(100. * dy_os / InfValue)
+    >>> from control import step_info, TransferFunction
+    >>> sys = TransferFunction([-1, 1], [1, 1, 1])
+    >>> S = step_info(sys)
+    >>> for k in S:
+    ...     print(f"{k}: {S[k]:3.4}")
+    ...
+    RiseTime: 1.256
+    SettlingTime: 9.071
+    SettlingMin: 0.9011
+    SettlingMax: 1.208
+    Overshoot: 20.85
+    Undershoot: 27.88
+    Peak: 1.208
+    PeakTime: 4.187
+    SteadyStateValue: 1.0
+
+    MIMO System: Simulate until a final time of 10. Get the step response
+    characteristics for the second input and specify a 5% error until the
+    signal is considered settled.
+
+    >>> from numpy import sqrt
+    >>> from control import step_info, StateSpace
+    >>> sys = StateSpace([[-1., -1.],
+    ...                   [1., 0.]],
+    ...                  [[-1./sqrt(2.), 1./sqrt(2.)],
+    ...                   [0, 0]],
+    ...                  [[sqrt(2.), -sqrt(2.)]],
+    ...                  [[0, 0]])
+    >>> S = step_info(sys, T=10., SettlingTimeThreshold=0.05)
+    >>> for k, v in S[0][1].items():
+    ...     print(f"{k}: {float(v):3.4}")
+    RiseTime: 1.212
+    SettlingTime: 6.061
+    SettlingMin: -1.209
+    SettlingMax: -0.9184
+    Overshoot: 20.87
+    Undershoot: 28.02
+    Peak: 1.209
+    PeakTime: 4.242
+    SteadyStateValue: -1.0
+    """
+    if isinstance(sysdata, (StateSpace, TransferFunction)):
+        if T is None or np.asarray(T).size == 1:
+            T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=True)
+        T, Yout = step_response(sysdata, T, squeeze=False)
+        if yfinal:
+            InfValues = np.atleast_2d(yfinal)
         else:
-            overshoot = 0
-
-        y_us = (np.sign(InfValue.real)*yout).min()
-        dy_us = np.abs(y_us)
-        if dy_us > 0:
-            undershoot = np.abs(100. * dy_us / InfValue)
+            InfValues = np.atleast_2d(sysdata.dcgain())
+        retsiso = sysdata.issiso()
+        noutputs = sysdata.noutputs
+        ninputs = sysdata.ninputs
+    else:
+        # Time series of response data
+        errmsg = ("`sys` must be a LTI system, or time response data"
+                  " with a shape following the python-control"
+                  " time series data convention.")
+        try:
+            Yout = np.array(sysdata, dtype=float)
+        except ValueError:
+            raise ValueError(errmsg)
+        if Yout.ndim == 1 or (Yout.ndim == 2 and Yout.shape[0] == 1):
+            Yout = Yout[np.newaxis, np.newaxis, :]
+            retsiso = True
+        elif Yout.ndim == 3:
+            retsiso = False
         else:
-            undershoot = 0
-
-        # RiseTime
-        rise_time = T[tr_upper_index] - T[tr_lower_index]
-
-        settling_max = (yout[tr_upper_index:]).max()
-        settling_min = (yout[tr_upper_index:]).min()
-
-    return {
-        'RiseTime': rise_time,
-        'SettlingTime': settling_time,
-        'SettlingMin': settling_min,
-        'SettlingMax': settling_max,
-        'Overshoot': overshoot,
-        'Undershoot': undershoot,
-        'Peak': peak_value,
-        'PeakTime': peak_time,
-        'SteadyStateValue': steady_state_value
-        }
+            raise ValueError(errmsg)
+        if T is None or Yout.shape[2] != len(np.squeeze(T)):
+            raise ValueError("For time response data, a matching time vector"
+                             " must be given")
+        T = np.squeeze(T)
+        noutputs = Yout.shape[0]
+        ninputs = Yout.shape[1]
+        InfValues = np.atleast_2d(yfinal) if yfinal else Yout[:, :, -1]
+
+    ret = []
+    for i in range(noutputs):
+        retrow = []
+        for j in range(ninputs):
+            yout = Yout[i, j, :]
+
+            # Steady state value
+            InfValue = InfValues[i, j]
+            sgnInf = np.sign(InfValue.real)
+
+            rise_time: float = np.NaN
+            settling_time: float = np.NaN
+            settling_min: float = np.NaN
+            settling_max: float = np.NaN
+            peak_value: float = np.Inf
+            peak_time: float = np.Inf
+            undershoot: float = np.NaN
+            overshoot: float = np.NaN
+            steady_state_value: complex = np.NaN
+
+            if not np.isnan(InfValue) and not np.isinf(InfValue):
+                # RiseTime
+                tr_lower_index = np.where(
+                    sgnInf * (yout - RiseTimeLimits[0] * InfValue) >= 0
+                    )[0][0]
+                tr_upper_index = np.where(
+                    sgnInf * (yout - RiseTimeLimits[1] * InfValue) >= 0
+                    )[0][0]
+                rise_time = T[tr_upper_index] - T[tr_lower_index]
+
+                # SettlingTime
+                settled = np.where(
+                    np.abs(yout/InfValue-1) >= SettlingTimeThreshold)[0][-1]+1
+                # MIMO systems can have unsettled channels without infinite
+                # InfValue
+                if settled < len(T):
+                    settling_time = T[settled]
+
+                settling_min = (yout[tr_upper_index:]).min()
+                settling_max = (yout[tr_upper_index:]).max()
+
+                # Overshoot
+                y_os = (sgnInf * yout).max()
+                dy_os = np.abs(y_os) - np.abs(InfValue)
+                if dy_os > 0:
+                    overshoot = np.abs(100. * dy_os / InfValue)
+                else:
+                    overshoot = 0
+
+                # Undershoot
+                y_us = (sgnInf * yout).min()
+                dy_us = np.abs(y_us)
+                if dy_us > 0:
+                    undershoot = np.abs(100. * dy_us / InfValue)
+                else:
+                    undershoot = 0
+
+                # Peak
+                peak_index = np.abs(yout).argmax()
+                peak_value = np.abs(yout[peak_index])
+                peak_time = T[peak_index]
+
+                # SteadyStateValue
+                steady_state_value = InfValue.real
+
+            retij = {
+                'RiseTime': rise_time,
+                'SettlingTime': settling_time,
+                'SettlingMin': settling_min,
+                'SettlingMax': settling_max,
+                'Overshoot': overshoot,
+                'Undershoot': undershoot,
+                'Peak': peak_value,
+                'PeakTime': peak_time,
+                'SteadyStateValue': steady_state_value
+                }
+            retrow.append(retij)
+
+        ret.append(retrow)
+
+    return ret[0][0] if retsiso else ret
 
 def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
                      transpose=False, return_x=False, squeeze=None):
@@ -1287,16 +1381,18 @@ def _default_time_vector(sys, N=None, tfinal=None, is_step=True):
         # only need to use default_tfinal if not given; N is ignored.
         if tfinal is None:
             # for discrete time, change from ideal_tfinal if N too large/small
-            N = int(np.clip(ideal_tfinal/sys.dt, N_min_dt, N_max))# [N_min, N_max]
-            tfinal = sys.dt * N
+            # [N_min, N_max]
+            N = int(np.clip(np.ceil(ideal_tfinal/sys.dt)+1, N_min_dt, N_max))
+            tfinal = sys.dt * (N-1)
         else:
-            N = int(tfinal/sys.dt)
-            tfinal = N * sys.dt # make tfinal an integer multiple of sys.dt
+            N = int(np.ceil(tfinal/sys.dt)) + 1
+            tfinal = sys.dt * (N-1) # make tfinal an integer multiple of sys.dt
     else:
         if tfinal is None:
             # for continuous time, simulate to ideal_tfinal but limit N
             tfinal = ideal_tfinal
         if N is None:
-            N = int(np.clip(tfinal/ideal_dt, N_min_ct, N_max)) # N<-[N_min, N_max]
+            # [N_min, N_max]
+            N = int(np.clip(np.ceil(tfinal/ideal_dt)+1, N_min_ct, N_max))
 
-    return np.linspace(0, tfinal, N, endpoint=False)
+    return np.linspace(0, tfinal, N, endpoint=True)


Test patch:
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index a1246103..afeb71ef 100644
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -336,41 +336,41 @@ def test_zmore_stability_margins(tsys_zmore):
 @pytest.mark.parametrize(
     'cnum, cden, dt,'
     'ref,'
-    'rtol, poly_is_inaccurate',
+    'rtol',
     [( # gh-465
       [2], [1, 3, 2, 0], 1e-2,
-      [ 2.955761, 32.398492,  0.429535,  1.403725,  0.749367,  0.923898],
-      1e-5, True),
+      [2.9558, 32.390, 0.43584, 1.4037, 0.74951, 0.97079],
+      2e-3), # the gradient of the function reduces numerical precision
      ( # 2/(s+1)**3
       [2], [1, 3, 3, 1], .1,
       [3.4927, 65.4212, 0.5763, 1.6283, 0.76625, 1.2019],
-      1e-4, True),
-     ( # gh-523 a
+      1e-4),
+     ( # gh-523
       [1.1 * 4 * np.pi**2], [1, 2 * 0.2 * 2 * np.pi,  4 * np.pi**2], .05,
       [2.3842, 18.161, 0.26953, 11.712, 8.7478, 9.1504],
-      1e-4, False),
-     ( # gh-523 b
-       # H1 = w1**2 / (z**2 + 2*zt*w1 * z + w1**2)
-       # H2 = w2**2 / (z**2 + 2*zt*w2 * z + w2**2)
-       # H = H1 * H2
-       # w1 = 1, w2 = 100, zt = 0.5
-      [5e4], [1., 101., 10101., 10100., 10000.], 1e-3,
-      [18.8766, 26.3564, 0.406841, 9.76358, 2.32933, 2.55986],
-      1e-5, True),
+      1e-4),
      ])
-@pytest.mark.filterwarnings("error")
-def test_stability_margins_discrete(cnum, cden, dt,
-                                    ref,
-                                    rtol, poly_is_inaccurate):
+def test_stability_margins_discrete(cnum, cden, dt, ref, rtol):
     """Test stability_margins with discrete TF input"""
     tf = TransferFunction(cnum, cden).sample(dt)
-    if poly_is_inaccurate:
-        with pytest.warns(UserWarning, match="numerical inaccuracy in 'poly'"):
-            out = stability_margins(tf)
-        # cover the explicit frd branch and make sure it yields the same
-        # results as the fallback mechanism
-        out_frd = stability_margins(tf, method='frd')
-        assert_allclose(out, out_frd)
-    else:
-        out = stability_margins(tf)
+    out = stability_margins(tf, method='poly')
     assert_allclose(out, ref, rtol=rtol)
+
+
+def test_stability_margins_methods():
+    # the following system gives slightly inaccurate result for DT systems
+    # because of numerical issues
+    omegan = 1
+    zeta = 0.5
+    resonance = TransferFunction(omegan**2, [1, 2*zeta*omegan, omegan**2])
+    omegan2 = 100
+    resonance2 = TransferFunction(omegan2**2, [1, 2*zeta*omegan2, omegan2**2])
+    sys = 5 * resonance * resonance2
+    sysd = sys.sample(0.001, 'zoh')
+    """Test stability_margins() function with different methods"""
+    out = stability_margins(sysd, method='best')
+    # confirm getting reasonable results using FRD method
+    assert_allclose(
+        (18.876634740386308, 26.356358386241055, 0.40684127995261044,
+         9.763585494645046, 2.3293357226374805, 2.55985695034263),
+        stability_margins(sysd, method='frd'), rtol=1e-5)
diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py
index 6df2493c..14e9692c 100644
--- a/control/tests/sisotool_test.py
+++ b/control/tests/sisotool_test.py
@@ -68,8 +68,8 @@ class TestSisotool:
 
         # Check the step response before moving the point
         step_response_original = np.array(
-            [0.    , 0.021 , 0.124 , 0.3146, 0.5653, 0.8385, 1.0969, 1.3095,
-             1.4549, 1.5231])
+            [0.    , 0.0216, 0.1271, 0.3215, 0.5762, 0.8522, 1.1114, 1.3221,
+             1.4633, 1.5254])
         assert_array_almost_equal(
             ax_step.lines[0].get_data()[1][:10], step_response_original, 4)
 
@@ -113,8 +113,8 @@ class TestSisotool:
 
         # Check if the step response has changed
         step_response_moved = np.array(
-            [0.    , 0.023 , 0.1554, 0.4401, 0.8646, 1.3722, 1.875 , 2.2709,
-              2.4633, 2.3827])
+            [0.    , 0.0237, 0.1596, 0.4511, 0.884 , 1.3985, 1.9031, 2.2922,
+             2.4676, 2.3606])
         assert_array_almost_equal(
             ax_step.lines[0].get_data()[1][:10], step_response_moved, 4)
 
@@ -157,7 +157,3 @@ class TestSisotool:
         # but 2 input, 1 output should
         with pytest.raises(ControlMIMONotImplemented):
             sisotool(sys221)
-
-
-
-
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 1eec5ead..983b9d7a 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -8,7 +8,6 @@ BG,  26 Jul 2020 merge statesp_array_test.py differences into statesp_test.py
 """
 
 import numpy as np
-from numpy.testing import assert_array_almost_equal
 import pytest
 import operator
 from numpy.linalg import solve
@@ -48,17 +47,6 @@ class TestStateSpace:
         """3-states square system (2 inputs x 2 outputs)"""
         return StateSpace(*sys322ABCD)
 
-    @pytest.fixture
-    def sys121(self):
-        """2 state, 1 input, 1 output (siso)  system"""
-        A121 = [[4., 1.],
-                [2., -3]]
-        B121 = [[5.],
-                [-3.]]
-        C121 = [[2., -4]]
-        D121 = [[3.]]
-        return StateSpace(A121, B121, C121, D121)
-
     @pytest.fixture
     def sys222(self):
         """2-states square system (2 inputs x 2 outputs)"""
@@ -763,70 +751,6 @@ class TestStateSpace:
             np.squeeze(sys322.horner(1.j)),
             mag[:, :, 0] * np.exp(1.j * phase[:, :, 0]))
 
-    @pytest.mark.parametrize('x',
-        [[1, 1], [[1], [1]], np.atleast_2d([1,1]).T])
-    @pytest.mark.parametrize('u', [0, 1, np.atleast_1d(2)])
-    def test_dynamics_and_output_siso(self, x, u, sys121):
-        assert_array_almost_equal(
-            sys121.dynamics(0, x, u),
-            sys121.A.dot(x).reshape((-1,)) + sys121.B.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys121.output(0, x, u),
-            sys121.C.dot(x).reshape((-1,)) + sys121.D.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys121.dynamics(0, x),
-            sys121.A.dot(x).reshape((-1,)))
-        assert_array_almost_equal(
-            sys121.output(0, x),
-            sys121.C.dot(x).reshape((-1,)))
-
-    # too few and too many states and inputs
-    @pytest.mark.parametrize('x', [0, 1, [], [1, 2, 3], np.atleast_1d(2)])
-    def test_error_x_dynamics_and_output_siso(self, x, sys121):
-        with pytest.raises(ValueError):
-            sys121.dynamics(0, x)
-        with pytest.raises(ValueError):
-            sys121.output(0, x)
-    @pytest.mark.parametrize('u', [[1, 1], np.atleast_1d((2, 2))])
-    def test_error_u_dynamics_output_siso(self, u, sys121):
-        with pytest.raises(ValueError):
-            sys121.dynamics(0, 1, u)
-        with pytest.raises(ValueError):
-            sys121.output(0, 1, u)
-
-    @pytest.mark.parametrize('x',
-        [[1, 1], [[1], [1]], np.atleast_2d([1,1]).T])
-    @pytest.mark.parametrize('u',
-        [[1, 1], [[1], [1]], np.atleast_2d([1,1]).T])
-    def test_dynamics_and_output_mimo(self, x, u, sys222):
-        assert_array_almost_equal(
-            sys222.dynamics(0, x, u),
-            sys222.A.dot(x).reshape((-1,)) + sys222.B.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys222.output(0, x, u),
-            sys222.C.dot(x).reshape((-1,)) + sys222.D.dot(u).reshape((-1,)))
-        assert_array_almost_equal(
-            sys222.dynamics(0, x),
-            sys222.A.dot(x).reshape((-1,)))
-        assert_array_almost_equal(
-            sys222.output(0, x),
-            sys222.C.dot(x).reshape((-1,)))
-
-    # too few and too many states and inputs
-    @pytest.mark.parametrize('x', [0, 1, [1, 1, 1]])
-    def test_error_x_dynamics_mimo(self, x, sys222):
-        with pytest.raises(ValueError):
-            sys222.dynamics(0, x)
-        with pytest.raises(ValueError):
-            sys222.output(0, x)
-    @pytest.mark.parametrize('u', [1, [1, 1, 1]])
-    def test_error_u_dynamics_mimo(self, u, sys222):
-        with pytest.raises(ValueError):
-            sys222.dynamics(0, (1, 1), u)
-        with pytest.raises(ValueError):
-            sys222.output(0, (1, 1), u)
-
-
 class TestRss:
     """These are tests for the proper functionality of statesp.rss."""
 
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 8705f380..a576d090 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -15,21 +15,20 @@ import numpy as np
 import pytest
 import scipy as sp
 
-
 import control as ct
-from control import (StateSpace, TransferFunction, c2d, isctime, isdtime,
-                     ss2tf, tf2ss)
+from control import StateSpace, TransferFunction, c2d, isctime, ss2tf, tf2ss
+from control.exception import slycot_check
+from control.tests.conftest import slycotonly
 from control.timeresp import (_default_time_vector, _ideal_tfinal_and_dt,
                               forced_response, impulse_response,
                               initial_response, step_info, step_response)
-from control.tests.conftest import slycotonly
-from control.exception import slycot_check
 
 
 class TSys:
     """Struct of test system"""
-    def __init__(self, sys=None):
+    def __init__(self, sys=None, call_kwargs=None):
         self.sys = sys
+        self.kwargs = call_kwargs if call_kwargs else {}
 
     def __repr__(self):
         """Show system when debugging"""
@@ -66,8 +65,8 @@ class TestTimeresp:
         T.initial = siso_ss1.yinitial - 9
         T.yimpulse = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766,
                                31.7344, 26.1668, 21.6292, 17.9245, 14.8945])
-        return T
 
+        return T
 
     @pytest.fixture
     def siso_tf1(self):
@@ -193,35 +192,150 @@ class TestTimeresp:
     def no_pole_cancellation(self):
         return TransferFunction([1.881e+06],
                                 [188.1, 1.881e+06])
-    
+
     @pytest.fixture
     def siso_tf_type1(self):
         # System Type 1 - Step response not stationary:  G(s)=1/s(s+1)
-        return TransferFunction(1, [1, 1, 0])
+        T = TSys(TransferFunction(1, [1, 1, 0]))
+        T.step_info = {
+             'RiseTime': np.NaN,
+             'SettlingTime': np.NaN,
+             'SettlingMin': np.NaN,
+             'SettlingMax': np.NaN,
+             'Overshoot': np.NaN,
+             'Undershoot': np.NaN,
+             'Peak': np.Inf,
+             'PeakTime': np.Inf,
+             'SteadyStateValue': np.NaN}
+        return T
 
     @pytest.fixture
     def siso_tf_kpos(self):
-        # SISO under shoot response and positive final value G(s)=(-s+1)/(s²+s+1)
-        return TransferFunction([-1, 1], [1, 1, 1])
+        # SISO under shoot response and positive final value
+        # G(s)=(-s+1)/(s²+s+1)
+        T = TSys(TransferFunction([-1, 1], [1, 1, 1]))
+        T.step_info = {
+             'RiseTime': 1.242,
+             'SettlingTime': 9.110,
+             'SettlingMin': 0.90,
+             'SettlingMax': 1.208,
+             'Overshoot': 20.840,
+             'Undershoot': 28.0,
+             'Peak': 1.208,
+             'PeakTime': 4.282,
+             'SteadyStateValue': 1.0}
+        return T
 
     @pytest.fixture
     def siso_tf_kneg(self):
-        # SISO under shoot response and negative final value k=-1 G(s)=-(-s+1)/(s²+s+1)
-        return TransferFunction([1, -1], [1, 1, 1])
+        # SISO under shoot response and negative final value
+        # k=-1 G(s)=-(-s+1)/(s²+s+1)
+        T = TSys(TransferFunction([1, -1], [1, 1, 1]))
+        T.step_info = {
+             'RiseTime': 1.242,
+             'SettlingTime': 9.110,
+             'SettlingMin': -1.208,
+             'SettlingMax': -0.90,
+             'Overshoot': 20.840,
+             'Undershoot': 28.0,
+             'Peak': 1.208,
+             'PeakTime': 4.282,
+             'SteadyStateValue': -1.0}
+        return T
 
     @pytest.fixture
-    def tf1_matlab_help(self):
-        # example from matlab online help https://www.mathworks.com/help/control/ref/stepinfo.html
-        return TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])
+    def siso_tf_step_matlab(self):
+        # example from matlab online help
+        # https://www.mathworks.com/help/control/ref/stepinfo.html
+        T = TSys(TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2]))
+        T.step_info = {
+            'RiseTime': 3.8456,
+            'SettlingTime': 27.9762,
+            'SettlingMin': 2.0689,
+            'SettlingMax': 2.6873,
+            'Overshoot': 7.4915,
+            'Undershoot': 0,
+            'Peak': 2.6873,
+            'PeakTime': 8.0530,
+            'SteadyStateValue': 2.5}
+        return T
+
+    @pytest.fixture
+    def mimo_ss_step_matlab(self):
+        A = [[0.68, -0.34],
+             [0.34, 0.68]]
+        B = [[0.18, -0.05],
+             [0.04, 0.11]]
+        C = [[0, -1.53],
+             [-1.12, -1.10]]
+        D = [[0, 0],
+             [0.06, -0.37]]
+        T = TSys(StateSpace(A, B, C, D, 0.2))
+        T.kwargs['step_info'] = {'T': 4.6}
+        T.step_info = [[{'RiseTime': 0.6000,
+                         'SettlingTime': 3.0000,
+                         'SettlingMin': -0.5999,
+                         'SettlingMax': -0.4689,
+                         'Overshoot': 15.5072,
+                         'Undershoot': 0.,
+                         'Peak': 0.5999,
+                         'PeakTime': 1.4000,
+                         'SteadyStateValue': -0.5193},
+                        {'RiseTime': 0.,
+                         'SettlingTime': 3.6000,
+                         'SettlingMin': -0.2797,
+                         'SettlingMax': -0.1043,
+                         'Overshoot': 118.9918,
+                         'Undershoot': 0,
+                         'Peak': 0.2797,
+                         'PeakTime': .6000,
+                         'SteadyStateValue': -0.1277}],
+                       [{'RiseTime': 0.4000,
+                         'SettlingTime': 2.8000,
+                         'SettlingMin': -0.6724,
+                         'SettlingMax': -0.5188,
+                         'Overshoot': 24.6476,
+                         'Undershoot': 11.1224,
+                         'Peak': 0.6724,
+                         'PeakTime': 1,
+                         'SteadyStateValue': -0.5394},
+                         {'RiseTime': 0.0000, # (*)
+                         'SettlingTime': 3.4000,
+                         'SettlingMin': -0.1034,
+                         'SettlingMax': -0.1485,
+                         'Overshoot': 132.0170,
+                         'Undershoot': 79.222, # 0. in MATLAB
+                         'Peak': 0.4350,
+                         'PeakTime': .2,
+                         'SteadyStateValue': -0.1875}]]
+                         # (*): MATLAB gives 0.4 here, but it is unclear what
+                         # 10% and 90% of the steady state response mean, when
+                         # the step for this channel does not start a 0 for
+                         # 0 initial conditions
+        return T
 
     @pytest.fixture
-    def tf2_matlab_help(self):
-        A = [[0.68, - 0.34], [0.34, 0.68]]
-        B = [[0.18], [0.04]]
-        C = [-1.12, - 1.10]
-        D = [0.06]
-        sys = StateSpace(A, B, C, D, 0.2)
-        return sys
+    def siso_ss_step_matlab(self, mimo_ss_step_matlab):
+        T = copy(mimo_ss_step_matlab)
+        T.sys = T.sys[1, 0]
+        T.step_info = T.step_info[1][0]
+        return T
+
+    @pytest.fixture
+    def mimo_tf_step_info(self,
+                          siso_tf_kpos, siso_tf_kneg,
+                          siso_tf_step_matlab):
+        Ta = [[siso_tf_kpos, siso_tf_kneg, siso_tf_step_matlab],
+               [siso_tf_step_matlab, siso_tf_kpos, siso_tf_kneg]]
+        T = TSys(TransferFunction(
+            [[Ti.sys.num[0][0] for Ti in Tr] for Tr in Ta],
+            [[Ti.sys.den[0][0] for Ti in Tr] for Tr in Ta]))
+        T.step_info = [[Ti.step_info for Ti in Tr] for Tr in Ta]
+        # enforce enough sample points for all channels (they have different
+        # characteristics)
+        T.kwargs['step_info'] = {'T_num': 2000}
+        return T
+
 
     @pytest.fixture
     def tsystem(self,
@@ -232,8 +346,9 @@ class TestTimeresp:
                 siso_dss1, siso_dss2,
                 mimo_dss1, mimo_dss2, mimo_dtf1,
                 pole_cancellation, no_pole_cancellation, siso_tf_type1,
-                siso_tf_kpos, siso_tf_kneg, tf1_matlab_help,
-                tf2_matlab_help):
+                siso_tf_kpos, siso_tf_kneg,
+                siso_tf_step_matlab, siso_ss_step_matlab,
+                mimo_ss_step_matlab, mimo_tf_step_info):
         systems = {"siso_ss1": siso_ss1,
                    "siso_ss2": siso_ss2,
                    "siso_tf1": siso_tf1,
@@ -254,8 +369,10 @@ class TestTimeresp:
                    "siso_tf_type1": siso_tf_type1,
                    "siso_tf_kpos": siso_tf_kpos,
                    "siso_tf_kneg": siso_tf_kneg,
-                   "tf1_matlab_help": tf1_matlab_help,
-                   "tf2_matlab_help": tf2_matlab_help,
+                   "siso_tf_step_matlab": siso_tf_step_matlab,
+                   "siso_ss_step_matlab": siso_ss_step_matlab,
+                   "mimo_ss_step_matlab": mimo_ss_step_matlab,
+                   "mimo_tf_step": mimo_tf_step_info,
                    }
         return systems[request.param]
 
@@ -309,102 +426,110 @@ class TestTimeresp:
         t, y = step_response(sys)
         np.testing.assert_array_equal(y, np.ones(len(t)))
 
-    def test_step_info(self):
-        # From matlab docs:
-        sys = TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])
-        Strue = {
-            'RiseTime': 3.8456,
-            'SettlingTime': 27.9762,
-            'SettlingMin': 2.0689,
-            'SettlingMax': 2.6873,
-            'Overshoot': 7.4915,
-            'Undershoot': 0,
-            'Peak': 2.6873,
-            'PeakTime': 8.0530,
-            'SteadyStateValue': 2.50
-        }
-
-        S = step_info(sys)
-
-        Sk = sorted(S.keys())
-        Sktrue = sorted(Strue.keys())
-        assert Sk == Sktrue
-        # Very arbitrary tolerance because I don't know if the
-        # response from the MATLAB is really that accurate.
-        # maybe it is a good idea to change the Strue to match
-        # but I didn't do it because I don't know if it is
-        # accurate either...
-        rtol = 2e-2
-        np.testing.assert_allclose([S[k] for k in Sk],
-                                   [Strue[k] for k in Sktrue],
-                                   rtol=rtol)
-
-    # tolerance for all parameters could be wrong for some systems
-    # discrete systems time parameters tolerance could be +/-dt
+    def assert_step_info_match(self, sys, info, info_ref):
+        """Assert reasonable step_info accuracy."""
+        if sys.isdtime(strict=True):
+            dt = sys.dt
+        else:
+            _, dt = _ideal_tfinal_and_dt(sys, is_step=True)
+
+        for k in ['RiseTime', 'SettlingTime', 'PeakTime']:
+            np.testing.assert_allclose(info[k], info_ref[k], atol=dt,
+                                       err_msg=f"{k} does not match")
+        for k in ['Overshoot', 'Undershoot', 'Peak', 'SteadyStateValue']:
+            np.testing.assert_allclose(info[k], info_ref[k], rtol=5e-3,
+                                       err_msg=f"{k} does not match")
+
+        # steep gradient right after RiseTime
+        absrefinf = np.abs(info_ref['SteadyStateValue'])
+        if info_ref['RiseTime'] > 0:
+            y_next_sample_max = 0.8*absrefinf/info_ref['RiseTime']*dt
+        else:
+            y_next_sample_max = 0
+        for k in ['SettlingMin', 'SettlingMax']:
+            if (np.abs(info_ref[k]) - 0.9 * absrefinf) > y_next_sample_max:
+                # local min/max peak well after signal has risen
+                np.testing.assert_allclose(info[k], info_ref[k], rtol=1e-3)
+
     @pytest.mark.parametrize(
-        "tsystem, info_true, tolerance",
-        [("tf1_matlab_help", {
-            'RiseTime': 3.8456,
-            'SettlingTime': 27.9762,
-            'SettlingMin': 2.0689,
-            'SettlingMax': 2.6873,
-            'Overshoot': 7.4915,
-            'Undershoot': 0,
-            'Peak': 2.6873,
-            'PeakTime': 8.0530,
-            'SteadyStateValue': 2.5}, 2e-2),
-         ("tf2_matlab_help", {
-             'RiseTime': 0.4000,
-             'SettlingTime': 2.8000,
-             'SettlingMin': -0.6724,
-             'SettlingMax': -0.5188,
-             'Overshoot': 24.6476,
-             'Undershoot': 11.1224,
-             'Peak': 0.6724,
-             'PeakTime': 1,
-             'SteadyStateValue': -0.5394}, .2),
-         ("siso_tf_kpos", {
-             'RiseTime': 1.242,
-             'SettlingTime': 9.110,
-             'SettlingMin': 0.950,
-             'SettlingMax': 1.208,
-             'Overshoot': 20.840,
-             'Undershoot': 27.840,
-             'Peak': 1.208,
-             'PeakTime': 4.282,
-             'SteadyStateValue': 1.0}, 2e-2),
-         ("siso_tf_kneg", {
-             'RiseTime': 1.242,
-             'SettlingTime': 9.110,
-             'SettlingMin': -1.208,
-             'SettlingMax': -0.950,
-             'Overshoot': 20.840,
-             'Undershoot': 27.840,
-             'Peak': 1.208,
-             'PeakTime': 4.282,
-             'SteadyStateValue': -1.0}, 2e-2),
-         ("siso_tf_type1", {'RiseTime': np.NaN,
-             'SettlingTime': np.NaN,
-             'SettlingMin': np.NaN,
-             'SettlingMax': np.NaN,
-             'Overshoot': np.NaN,
-             'Undershoot': np.NaN,
-             'Peak': np.Inf,
-             'PeakTime': np.Inf,
-             'SteadyStateValue': np.NaN}, 2e-2)],
+        "yfinal", [True, False], ids=["yfinal", "no yfinal"])
+    @pytest.mark.parametrize(
+        "systype, time_2d",
+        [("ltisys", False),
+         ("time response", False),
+         ("time response", True),
+         ],
+        ids=["ltisys", "time response (n,)", "time response (1,n)"])
+    @pytest.mark.parametrize(
+        "tsystem",
+        ["siso_tf_step_matlab",
+         "siso_ss_step_matlab",
+         "siso_tf_kpos",
+         "siso_tf_kneg",
+         "siso_tf_type1"],
         indirect=["tsystem"])
-    def test_step_info(self, tsystem, info_true, tolerance):
-        """  """
-        info = step_info(tsystem)
+    def test_step_info(self, tsystem, systype, time_2d, yfinal):
+        """Test step info for SISO systems."""
+        step_info_kwargs = tsystem.kwargs.get('step_info', {})
+        if systype == "time response":
+            # simulate long enough for steady state value
+            tfinal = 3 * tsystem.step_info['SettlingTime']
+            if np.isnan(tfinal):
+                pytest.skip("test system does not settle")
+            t, y = step_response(tsystem.sys, T=tfinal, T_num=5000)
+            sysdata = y
+            step_info_kwargs['T'] = t[np.newaxis, :] if time_2d else t
+        else:
+            sysdata = tsystem.sys
+        if yfinal:
+            step_info_kwargs['yfinal'] = tsystem.step_info['SteadyStateValue']
 
-        info_true_sorted = sorted(info_true.keys())
-        info_sorted = sorted(info.keys())
+        info = step_info(sysdata, **step_info_kwargs)
 
-        assert info_sorted == info_true_sorted
+        self.assert_step_info_match(tsystem.sys, info, tsystem.step_info)
 
-        np.testing.assert_allclose([info_true[k] for k in info_true_sorted],
-                                   [info[k] for k in info_sorted],
-                                   rtol=tolerance)
+    @pytest.mark.parametrize(
+        "yfinal", [True, False], ids=["yfinal", "no_yfinal"])
+    @pytest.mark.parametrize(
+        "systype", ["ltisys", "time response"])
+    @pytest.mark.parametrize(
+        "tsystem",
+        ['mimo_ss_step_matlab',
+         pytest.param('mimo_tf_step', marks=slycotonly)],
+        indirect=["tsystem"])
+    def test_step_info_mimo(self, tsystem, systype, yfinal):
+        """Test step info for MIMO systems."""
+        step_info_kwargs = tsystem.kwargs.get('step_info', {})
+        if systype == "time response":
+            tfinal = 3 * max([S['SettlingTime']
+                              for Srow in tsystem.step_info for S in Srow])
+            t, y = step_response(tsystem.sys, T=tfinal, T_num=5000)
+            sysdata = y
+            step_info_kwargs['T'] = t
+        else:
+            sysdata = tsystem.sys
+        if yfinal:
+            step_info_kwargs['yfinal'] = [[S['SteadyStateValue']
+                                           for S in Srow]
+                                          for Srow in tsystem.step_info]
+
+        info_dict = step_info(sysdata, **step_info_kwargs)
+
+        for i, row in enumerate(info_dict):
+            for j, info in enumerate(row):
+                self.assert_step_info_match(tsystem.sys,
+                                            info, tsystem.step_info[i][j])
+
+    def test_step_info_invalid(self):
+        """Call step_info with invalid parameters."""
+        with pytest.raises(ValueError, match="time series data convention"):
+            step_info(["not numeric data"])
+        with pytest.raises(ValueError, match="time series data convention"):
+            step_info(np.ones((10, 15)))                     # invalid shape
+        with pytest.raises(ValueError, match="matching time vector"):
+            step_info(np.ones(15), T=np.linspace(0, 1, 20))  # time too long
+        with pytest.raises(ValueError, match="matching time vector"):
+            step_info(np.ones((2, 2, 15)))                   # no time vector
 
     def test_step_pole_cancellation(self, pole_cancellation,
                                     no_pole_cancellation):
@@ -412,13 +537,9 @@ class TestTimeresp:
         # https://github.com/python-control/python-control/issues/440
         step_info_no_cancellation = step_info(no_pole_cancellation)
         step_info_cancellation = step_info(pole_cancellation)
-        for key in step_info_no_cancellation:
-            if key == 'Overshoot':
-                # skip this test because these systems have no overshoot
-                # =>  very sensitive to parameters
-                continue
-            np.testing.assert_allclose(step_info_no_cancellation[key],
-                                       step_info_cancellation[key], rtol=1e-4)
+        self.assert_step_info_match(no_pole_cancellation,
+                                    step_info_no_cancellation,
+                                    step_info_cancellation)
 
     @pytest.mark.parametrize(
         "tsystem, kwargs",
@@ -634,20 +755,23 @@ class TestTimeresp:
          (TransferFunction(1, [1, .5, 0]), 25)])         # poles at 0.5 and 0
     def test_auto_generated_time_vector_tfinal(self, tfsys, tfinal):
         """Confirm a TF with a pole at p simulates for tfinal seconds"""
-        np.testing.assert_almost_equal(
-            _ideal_tfinal_and_dt(tfsys)[0], tfinal, decimal=4)
+        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(tfsys)
+        np.testing.assert_allclose(ideal_tfinal, tfinal, rtol=1e-4)
+        T = _default_time_vector(tfsys)
+        np.testing.assert_allclose(T[-1], tfinal, atol=0.5*ideal_dt)
 
     @pytest.mark.parametrize("wn, zeta", [(10, 0), (100, 0), (100, .1)])
-    def test_auto_generated_time_vector_dt_cont(self, wn, zeta):
+    def test_auto_generated_time_vector_dt_cont1(self, wn, zeta):
         """Confirm a TF with a natural frequency of wn rad/s gets a
         dt of 1/(ratio*wn)"""
 
         dtref = 0.25133 / wn
 
         tfsys = TransferFunction(1, [1, 2*zeta*wn, wn**2])
-        np.testing.assert_almost_equal(_ideal_tfinal_and_dt(tfsys)[1], dtref)
+        np.testing.assert_almost_equal(_ideal_tfinal_and_dt(tfsys)[1], dtref,
+                                       decimal=5)
 
-    def test_auto_generated_time_vector_dt_cont(self):
+    def test_auto_generated_time_vector_dt_cont2(self):
         """A sampled tf keeps its dt"""
         wn = 100
         zeta = .1
@@ -674,21 +798,23 @@ class TestTimeresp:
     def test_default_timevector_functions_c(self, fun):
         """Test that functions can calculate the time vector automatically"""
         sys = TransferFunction(1, [1, .5, 0])
+        _tfinal, _dt = _ideal_tfinal_and_dt(sys)
 
         # test impose number of time steps
         tout, _ = fun(sys, T_num=10)
         assert len(tout) == 10
 
         # test impose final time
-        tout, _ = fun(sys, 100)
-        np.testing.assert_allclose(tout[-1], 100., atol=0.5)
+        tout, _ = fun(sys, T=100.)
+        np.testing.assert_allclose(tout[-1], 100., atol=0.5*_dt)
 
     @pytest.mark.parametrize("fun", [step_response,
                                      impulse_response,
                                      initial_response])
-    def test_default_timevector_functions_d(self, fun):
+    @pytest.mark.parametrize("dt", [0.1, 0.112])
+    def test_default_timevector_functions_d(self, fun, dt):
         """Test that functions can calculate the time vector automatically"""
-        sys = TransferFunction(1, [1, .5, 0], 0.1)
+        sys = TransferFunction(1, [1, .5, 0], dt)
 
         # test impose number of time steps is ignored with dt given
         tout, _ = fun(sys, T_num=15)
@@ -696,7 +822,7 @@ class TestTimeresp:
 
         # test impose final time
         tout, _ = fun(sys, 100)
-        np.testing.assert_allclose(tout[-1], 100., atol=0.5)
+        np.testing.assert_allclose(tout[-1], 100., atol=0.5*dt)
 
 
     @pytest.mark.parametrize("tsystem",