I am modeling an analog filter with digital software and have reduced the model to a 4th order FIR filter in discrete space with transfer function
$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + b_4 z^{-4}}{1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + a_4 z^{-4}} $$
It happens to be unstable, but I have noticed that most analog "musical" LP filters distort non-linearly when in this realm of high resonance rather than blow up to infinity (of course). Is it possible to simulate this by adding a nonlinear element to the IIR block diagram?
Clipping with $y[n] = \operatorname{clip}(x[n] + \ldots)$ creates a periodic [-1, -1, 1, 1, ...] in the output signal, so is there another way I can prevent numerical blowup?