Reconstructed output mismatch for LTI system
The problem is that the formula $$H(f) = \frac{Y(f)}{X(f)},$$ where $H(f)$ is the frequency response, and $X(f)$, $Y(f)$ are the input and output, is valid only for frquencies $f$ where $|X(f)| \neq 0$.

Furthermore, if $|X(f_0)| \approx 0$, you're going to run into numerical errors, since $|H(f_0)| \rightarrow \infty$.

Your input signal has a bunch of places where it's zero, or close to zero. I wrote a little function that takes care of those cases:

function h = impres(x,y)

X = fft(x);
Y = fft(y);

for i = 1:length(X)
    if abs(X(i)) &gt; 10
        H(i) = Y(i)/X(i);
    else
        H(i) = 0;
    end
end

h = ifft(H);

end


While not perfect, and could be tweaked, I get:



where the red line is the estimated output.

One more thing to be aware of: the equivalent in time of the IDFT is the circular convolution, not the linear convolution. In this case there is not too much difference, but you want to be as correct as you can in your code.