Fixes python-control#48

The old algorithm called a general-purpose ODE integrator at each
timestep, and this could be extremely slow (hundreds of times slower
than scipy.signal.lsim).  This new algorithm evaluates the matrix
exponential once, and uses it to advance the simulation at each step.

Importantly, the old algorithm would work for variable timesteps, but
the new version assumes a fixed timestep (as does Matlab's lsim).  This
is the usual case, and it is the price one pays for increased speed.

Note that this algorithm is the same idea as the algorithm used in
scipy.signal.lsim, but the scipy version requires inverting the `A`
matrix, and thus does not work if the system has a pole at the origin.
The algorithm used here works even if there is a pole at the origin.
This commit also adds a test for this case (the double integrator).