double[] eigenvalues = m.eig().getRealEigenvalues();
double[] eigRange = range(eigenvalues, minTimeConstant, maxTimeConstant);
//We will give an approximation of an impulse, then let the system decay for 10 X the longest time constant,
//then step in increments of 1/2 X the longest time constant until there is little change in the integral
Integrator integrator = new EulerIntegrator(-1f / (10f * (float) eigRange[0]));
float[] integral = new float[system.getOutputDimension()];
system.setState(new float[eigenvalues.length]);
float[] impulse = new float[system.getInputDimension()];
float pulseWidth = (float) minTimeConstant; //prevents saturation
float pulseHeight = 1f / pulseWidth;
impulse[input] = pulseHeight;
float[] zero = new float[system.getInputDimension()];
Units[] units = new Units[system.getInputDimension()];
TimeSeries pulse = integrator.integrate(system,
new TimeSeriesImpl(new float[]{0f, pulseWidth}, new float[][]{impulse, impulse}, units));
integral = integrate(pulse);
float decayTime = -10f / (float) eigRange[1];
TimeSeries decay = integrator.integrate(system,
new TimeSeriesImpl(new float[]{0f, decayTime}, new float[][]{zero, zero}, units)); //time-invariant, so we can start at 0
float[] increment = integrate(decay);
integral = MU.sum(integral, increment);
float stepTime = -.5f / (float) eigRange[1];
do {
decay = integrator.integrate(system,
new TimeSeriesImpl(new float[]{0f, stepTime}, new float[][]{zero, zero}, units));
increment = integrate(decay);
integral = MU.sum(integral, increment);
} while (substantialChange(integral, increment, .001f * stepTime / decayTime));