As pointed out in issue #239, there is an inconsistency between the way that continuous time simulations and discrete time simulations return their outputs for single output systems. In continuous time, the code follows the documentation and the following works:

t, y = step_response(sys)
plot(t, y)

However, in discrete time simulations the returned value for y had shape (1, len(t)), which meant that you had to do the following for discrete time systems:

t, y = step_response(sysd)
plot(t, y[0])

The problem was that for continuous time simulations forced_response() (the ultimate function that gets called) had an np.squeeze() at the end but this wasn't present for discrete time simulations. This PR fixes that problem so that the behavior is consistent, resolving the bug in issue #239.

Other small changes:

• Added unit tests that catch the original bug

• Added a squeeze keyword that controls whether this behavior is implemented (default = True).

• Did some PEP8 cleanup on timeresponse.py

Note that this may make some existing code break since discrete time simulations can now return an output that is a 1D array instead of 2D. The fix is either to remove the output index (y[0] becomes y) or set the squeeze keyword to False).