Pre-emphasis / De-emphasis filter design
You really can use any filter design scheme that you want. I can't think of any reason why you would need to use an IIR approach, exclusively, so you have a number of options. For FIR filters, you have choices like the frequency sampling and least squares methods (look at fir2 and firls in MATLAB). 

For an IIR implementation, you should be able to use the bilinear transform straightforwardly; since the filters you referenced are simple single-pole and single-zero ems, it's probably just as easy to work the math out by hand. One thing to remember is that you'll want to prewarp the frequency axis first to ensure that you get a good approximation across the passband of interest. It might help for you to provide any error that you were getting from MATLAB; the lack of poles shouldn't be problem.

Edit: A few things to address your comments:


I wouldn't say that firls, being a least-squares-based approach, is a step-wise approximation of a smooth curve. Yes, you provide a number of samples of the desired response to the least-squares algorithm, but you should still get a smooth curve out.
I started to work out the math on applying the bilinear transform to the pre-emphasis filter and came across a problem. Its frequenecy response is: 
$$
H_{pe}(f) = 1 + j2\pi f \tau_x
$$
Noting that for the frequency response $H_{pe}(f)$, $s = j\omega = j 2 \pi f$:
$$
H_{pe}(s) = 1 + s\tau_x
$$
The bilinear transform maps the $j\omega$ axis in the $s$-plane to the unit circle of the discrete-time approximation's $z$-plane. So, we should be concerned about what happens to the transfer function as $s$ is allowed to vary along the entire imaginary axis. Note that:
$$
\lim_{s \to \infty} H_{pe}(s) = \infty
$$
This means that, for the discrete-time approximation generated by the bilinear transform, as the angular frequency $\omega \to \pm \pi$, the discrete-time filter's magnitude response will blow up to $\infty$, which is probably not what you want.

This problem would present itself for any analog system whose transfer function has a larger-degree numerator than denominator. So it's probably justified why MATLAB issues the following error:

&gt;&gt; [num, den] = bilinear([75e-6 1], 1, 100e3);
Error using bilinear (line 69)
Numerator cannot be higher order than denominator.


(Note that in your example, you passed an empty vector for the denominator; instead, you should pass 1, since the denominator of the transfer function is unity)


So as it turns out, I'm not aware of any analog-mapping techniques that will definitely work for your application. I don't believe the impulse invariant method will work since the pre-emphasis filter isn't bandlimited. You could use one of the FIR design methods, an IIR least-squares technique (e.g. Berchin's frequency-domain least-squares (FDLS) algorithm), or perhaps look into the matched Z-transform method, which I'm not familiar with.