This PR allows parameter values to be saved in time responses as well as passing parameter values (via the params keyword) for the initial_response, ,step_response, and final_response functions (which call through to input_output_response for nonlinear I/O systems).

With this functionality, you can do things like

    nlsys = ct.nlsys(
        lambda t, x, u, params: params['a'] * x[0] + u[0],
        states = 1, inputs = 1, outputs = 1, params={'a': 0})
    timevec = np.linspace(0, 1)
    resp = step_response(nlsys, timevec, params={'a': -0.5})

The simulation will run with params['a'] = -0.5 and this information will be stored in resp.params.