The iirfilter() function internally generates filter prototypes in zpk format, transforms them to tf format, and then back to zpk format:

z, p, k = typefunc(N)
b, a = zpk2tf(z, p, k)
b, a = lp2lp(b, a, wo=warped)
b, a = bilinear(b, a, fs=fs)
return tf2zpk(b, a)

But conversion to tf format introduces numerical errors, causing higher-order filters to fail, even though the higher-order prototypes work fine. It should use zpk format (or state-space format?) throughout, and these functions should be changed to accept zpk format or be replaced by ones that do.

I started to implement this, but then realized it involves lots of changes to lots of functions, some of which I'm not sure the best way to do, so registering this issue in case others want to work on it, too.

For instance, matlab's bilinear function accepts tf, zpk, and ss formats, by varying the number of input parameters, which I don't think is "Pythonic" (though lti() does it). SciPy's bilinear function accepts b and a as separate parameters, and then fs, so it can't be modified to use variable number of input variables anyway, without breaking compatibility. A new function could be created that takes zpk, or the existing function could take a list as its first parameter instead, which could be zpk, tf, or ss, or something else.

For comparison, Octave uses sftrans instead of lp2lp, and only accepts zpk format.