Can you show the Matlab code you used to get the result you wanted? An ideal differentiator s is not proper tf so I do not know how c2d was used on it to covert it to discrete tf.

A tf should be ratio of 2 polynomials, hopefully with number of zeros less than or equal the number of poles. G(s)=s has more zeros than poles.

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.