The lsim command can be much, much slower than scipy.signal.lsim. For instance, consider the following example (run in IPython for %timeit):

import scipy.signal as sig
import control
import numpy as np
t = np.arange(0,1000,0.5)
u = np.sin(2.*t)
sys1 = sig.lti(1.,[1,0,1])
sys2 = control.tf(1.,[1,0,1])
%timeit sig.lsim(sys1,u,t)
%timeit control.lsim(sys2,u,t)

For me, scipy.signal.lsim gives 17.4 ms per loop, and control.lsim gives 7.64 s per loop.

The current implementation of control.lsim calls a general-purpose ODE solver at each step. I think it is better to assume a linear system with constant-spaced timesteps in t (as done in scipy.signal.lsim and in Matlab's lsim), and then one can evaluate the matrix exponential once and apply it at each step. If non-uniform timesteps are desired (which is very rare, at least in my experience), then we could have an alternative routine control.lsim2 that calls the general-purpose ODE solver--that is what is done in scipy.signal.