Fourier transform vs lsim, results dont match.

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

Y(e^jω)=G(e^jω)U(e^jω).

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

Implementation:

-U(e^jω) is obtained easily enough using fft(u).

-For G(e^jω) I use freqresp(sys,ω) (and build ω 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 yprime(t)= lsim(G,u). Unfortunately these dont match. What am I missing? The lsim results appear to have some integrator effect in them???

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')

[Niek]

Changed system G to a (stable) closed loop system, results seem to match a lot better now! Only some transient behaviour on the start of the lsim data