In most circumstances the bilinear transform is superior. The principal disadvantage of the BLT is that the frequency scale is ``warped'', i.e. distorted in a non-linear fashion, so that, e.g. for a LP filter,
|0 Hz||in the analogue domain maps on to||0 Hz||in the digital domain|
|fc||in the analogue domain maps on to||fc||in the digital domain|
|infinity Hz||in the analogue domain maps on to||fs/2||in the digital domain|
In most cases this is just what you want: for a LP or BP filter you want a response which is exactly zero at the Nyquist frequency (fs/2), and for a HP filter you want the response to be exactly unity at this frequency.
The one exception to this is the Bessel filter. The only advantage of a Bessel design over a Butterworth or Chebyshev is that the phase response is nearly linear throughout the passband, i.e. the group delay is almost constant throughout the passband. The ``warping'' inherent in the bilinear transform method upsets this linearity. To get round this problem, tick (check) the ``matched z-transform'' box. The matched z-transform does not warp the frequency scale, so a digital Bessel filter designed by this method will have the near-linear phase characteristic you'd like to see.
In fact, a low-pass filter designed by MZT is impulse-invariant with the analogue prototype, which means that it has the same impulse response, so it behaves identically to the corresponding analogue filter. Strictly speaking this is true only if you ignore aliasing: if the response of the digital LP filter is not negligible at fs/2, the response in the vicinity of this frequency will be perturbed by the spurious response at frequencies above fs/2, which are folded down into the band below fs/2. That's why you're advised to check the frequency response graphs carefully if you use this method!
A further advantage of the MZT is that the number of zeros is reduced, compared with the BLT case. Further information is given in the PostScript documentation (here is pdf version) for the mkfilter program. An LP filter designed by MZT has no z-plane zeros. This means that no xvals vector is required; this improves efficiency. Have a look at the generated code.
To summarize: you might use the matched z-transform if