function f = second_order_system(k,zeta,wn)

%simulate the response of the canonical second order LTI system
% 0 < zeta < 1
close all

if zeta > 0
    t = linspace(0,6/(zeta*wn),5000);
    T = 1/(zeta * wn);
else
    t = linspace(0,6*pi/wn,5000);
    T = NaN;
end

num = [k * wn^2];
den = [1, 2 * zeta * wn, wn^2];

y = step(num,den,t);
%wn = wn*sqrt(1-zeta^2);
el = k * (1 - exp(-zeta * wn * t));
eu = k * (1 + exp(-zeta * wn * t));
plot([-1 0 t],[0 0 y'],'LineWidth',2)
hold on
plot([-1 0 t],[0 0 ones(1,length(t))],'r','LineWidth',2)
plot(t,el,'k--','LineWidth',1);
plot(t,eu,'k--','LineWidth',1);
plot([0 T],[0 k],'k--','LineWidth',1)

set(gca,'Xlim',[-1 max(t)],'Ylim',sort([0 max(eu)]),'YTick',sort([0 1-eps k]), ...
    'XTick',[0 T],'XTicklabel',{'0','T'})
set(gca,'FontSize',16,'FontWeight','b','FontAngle','i')
grid on

xlabel('Time [s]')
th1 = text(0.8*max(t),1.25,'u(t)=h(t)','Color','r');
th2 = text(T,1.1*max(y),'y(t)','Color','b');
th3 = text(0.5*max(t),0.5*k,['\zeta = ' num2str(zeta)],'Color','k');
th4 = text(0.5*max(t),0.3*k,['\omega_n = ' num2str(wn)],'Color','k');
th5 = text(0.5*max(t),0.1*k,['k = ' num2str(k)],'Color','k');
set([th1 th2 th3 th4 th5],'FontSize',16,'FontWeight','b','FontAngle','i')