num [T_d*k k]
den [t^2 2*zeta*t 1]

G = tf (num, den)
t_d = 20;
G_d=c2d (G, t_d, 'zoh')
[num_d, den_d]=tfdata (G_d, 'v')

y_p1=0;
y_p2=0;
u_p1=0;
u_p2=0;

y_d=zeros(N,1)
t_d=zeros(N,1)
N=21;
u=2;
for i=1:N
	y=-den_d(2)*y_p1-den_d(3)*y_p2+num_d(2)*u_p1+num_d(3)*u_p2;
	y_p2=y_p1;
	y_p1=y;
	u_p2=u_p1;
	u_p1=u;
y_d(i)=y;
t_d(i)=(i-1)*T_d;
	end

[y_k, t_k]=step (G(N-1)*T_d);
plot (t_k, y_k, 'b');
hold on
grid on
stairs (t_d, y_d, 'r');
hold off