I just did a direct 1:1 translation between Matlab and Mathematica. No problem any where. For the final point you had with
sysz1=tf(numz,denz,ts,'variable','z^-1');
This is just a way to write the z transform in DSP convention. The tf does not change! All what Matlab does is divide by z^2 both numerator and denominator. You can see that it is the same tf like this:
sysz2=tf(numz,denz,ts);
sysz2DSP=tf(numz,denz,ts,'Variable','z^-1');
close all
[Y1,T1] = step(sysz2DSP);
[Y,T] = step(sysz2);
subplot(1,2,1); plot(T1,Y1); subplot(1,2,2); plot(T,Y)
You'll see it is the same response. Since it is the same transfer function. In Mathematica, here is the equivalent code to your Matlab's:
k = 10;
ksi1 = 0.01;
ksi2 = k*ksi1;
w1 = 200*2*Pi;
w2 = w1;
num = (s^2 + 2*ksi1*w1*s + w1^2);
den = (s^2 + 2*ksi2*w2*s + w1^2);
sys = TransferFunctionModel[ num/den, s]
ts = 0.000442;
numz = z^2 (1 - ts*ksi1*w1 + w1^2*ts^2/4 ) - (2 + w1^2*ts^2/2)*z + (1 + ts*ksi1*w1 + w1^2*ts^2/4);
denz = z^2 (1 - ts*ksi2*w2 + w2^2*ts^2/4) - z*(2 + w2^2*ts^2/2) + ( 1 + ts*ksi2*w2 + w2^2*ts^2/4);
sysz = TransferFunctionModel[ numz/denz, z, SamplingPeriod -> ts]
ts = 0.000442;
numz = z^2 (1 - ts*ksi1*w1 + w1^2*ts^2/4) - z*(2 + w1^2*ts^2/2) + ( 1 + ts*ksi1*w1 + w1^2*ts^2/4);
denz = z^2 (1 - ts*ksi2*w2 + w2^2*ts^2/4) - z*(2 + w2^2*ts^2/2 ) + (1 + ts*ksi2*w2 + w2^2*ts^2/4);
sysz2 = TransferFunctionModel[ numz/denz, z, SamplingPeriod -> ts]
Gives the same exact result as Matlab's. Again, What you did in Matlab using 1/z instead of z, does not change the transfer function. It is just another way of writing it.
