This PR implements MIMO impulse and step responses, mirroring the style and functionality of frequency response functions for MIMO systems. It resolves issue #512 and is consistent with the discussion issue #453. Specifically:

• If a system is MIMO, calling step_response or impulse_response will generate an array of responses (indexed by output, input, and time).
• Consistent with frequency response methods, if a system is SISO then by default the output will just be indexed by time.
• Input and output axes that are single dimensional can be removed by using the squeeze=True keyword.
• SISO systems can be tried as MIMO systems (retaining the output and input axes in the response) using the squeeze=False keyword.

Note that the response for initial_response and forced_response is not indexed by the input (since there is no input in one case and a specifically applied input in the other). The processing of the squeeze keyword is the same: by default SISO systems will be returned without the output axis, use squeeze=True to remove single dimensional axes from MIMO systems, use squeeze=False to force SISO systems to be treated like MIMO systems.

Other (small changes):

• Updated unit tests to make sure all of the squeezing rules are properly implemented with the new MIMO functionality.
• Found a bug in the way that the transpose keyword was working: it changed [output, input, time] to [time, input, output] instead of [time, output, input]. This was caught when adding unit tests.
• Found some errors in the MATLAB unit tests where the return arguments were in the wrong order (didn't show up with SISO only responses).