MATLAB: Fourier transform vs $\tt lsim$, results don't match My (simplified) problem:
Suppose I have a signal $u$, which passes through system/filter $G$, which gives output signal $y$. Then I have a theorem that states:

If a signal $u(t)$, with Fourier transform $U(e^{j \omega})$ is applied to system $G$, then the Fourier transform of the output $y(t)$ is given by:

$$Y(e^{j \omega})=G(e^{j \omega})U(e^{j \omega})$$

I try to verify this in MATLAB, but something goes wrong. 

Implementation: 


$U(e^{j \omega})$ is obtained easily enough using fft(u).
For $G(e^{j \omega})$ I use freqresp(sys,$\omega$) (and build
$\omega$ in such a way that it is in the same format as the fft
results: [dc positive-freqs negative-freqs]). I verify this by comparing the result to the bode plot (see commented lines in script)


To verify I perform the inverse Fourier transform y = ifft(Y) and compare the obtained $y(t)$ to the lsim result, which is obtained with $y_{prime}(t)$ = lsim(G,u). Unfortunately these dont match. What am I missing?

clear 
close all

G = tf([4700 4393 3.245e08],[1 7.574 1.202e5 0 0]);

f_sh = 1000;                                    %sampling frequency

Gd = c2d(G,1/f_sh);                             %plant discretization

t_end = 10;                                     %[s] simulation time end
t = linspace(0,t_end,t_end*f_sh+1);               %time vector
f_ref = 0.5;                                    %reference freq
u = sin(2*pi*f_ref*t);                        %reference signal

L = length(u);
dF = f_sh/L;                                    %frequency bin step

w = dF:dF:f_sh/2;                               %one-sided frequency vector in Hz

w_eval = 2*pi*[0 w -fliplr(w)];                 %frequency vector for freqresp, formatted to match fft results

H = freqresp(Gd,w_eval);                        %frequency response of system wrt freq vec

H = transpose(squeeze(H));                      %remove singelton dimensions

% bode(Gd) %should be equal to the semilog plot
% figure
% semilogx(w,20*log10(abs(H(2:(L+1)/2))))             %plot in Hz   


U = fft(u);
Y = H.*U;
y = ifft(Y);                                    %frequency domain response

yprime = lsim(Gd,u);                              %simulation/time domain resonse

figure
plot(t,yprime,t,y)                           %y should equal yprime
legend('lsim response','fft(sys)*fft(ref) response')


The lsim results appear to have some integrator effect in them???
To perform the FFT approach correctly, you will need the FFT length to be at least $N_U + N_H - 1$ where $N_U$ is the length of $u$ and $N_H$ is the length of $h$.  

It looks like $h$ is an IIR filter, so the FFT approach probably won't work well.

The Mathworks has a nice example comparing circular with linear convolution.

IIR filters are best implemented in the time domain. For their specs, they are usually much lower order than FIR filters and so we can afford to implement them directly.