Compared results of
sysc = ct.TransferFunction([1, z], [1, p])
sysd = ct.matlab.c2d(sysc, Ts, method='bilinear')
sysd1 = ct.matlab.c2d(sysc, Ts, method='matched')

yields a sysd1 with a gain that is about 10x that of sysd - is that expected?