Tend=1000;

fs=2100;                                   % sampling rate
t = 0:1/fs:Tend;                           % sampling instants
N=length(t);
 
input=[1;zeros(N-1,1)];                    % modelling the obtaning impulse reponse
ws=tf([10],[1 1 1 1 1]);
data=lsim(ws,input,t);
 
temp = fft(data,N);                        % perform N-point transform
fftdata = temp(1:N/2);                     % just look at first half                              

 
m = 20*log10(abs(fftdata));                % m = magnitude of sinusoids
p = unwrap(angle(fftdata));                % p = phase of sinusoids, unwrap()
                                           % copes with 360 degree jumps

f = (0:N/2-1)*fs/(N);                      % calculate Hertz values

% plot spectrum 0..fs/2 :
figure('Name','MyBode')
subplot(2,1,1), semilogx(f,m);             % plot magnitudes
ylabel('Abs. Magnitude'), grid on;
subplot(2,1,2), semilogx(f,p*180/pi);      % plot phase in degrees
ylabel('Phase [Degrees]'), grid on;
xlabel('Frequency [Hertz]');

figure('Name','Bode'),bode(ws);
 
clear