Force an unstable IIR filter to be stable by forcing a nonlinearity in the digital block diagram
You won't be able to avoid the instability with a linear model. With your current discrete model, the filter will blow up at a certain resonance setting. Using a better filter topology or higher precision will allow you to go a bit farther, but it will blow up.

Preventing that would require a saturating nonlinearity in the feedback path. That path is hidden in your implementation. To access it, you need to go back to the 4-section model without feedback and add the nonlinearity when you add the feedback path. But then you won't be able to simplify the filter to a fourth order linear filter.

Your options are then to delay the feedback path with its nonlinearity by one sample and feed it back into the 4-section model input. That will create significant deviations from the intended filter behaviour for high frequencies and high resonances. You can compensate for that somewhat, but it is tedious to get right. See Antti Huovilainen's paper for an approach like that.

The other option is to solve the implicit feedback equation using an iterative method for every single sample. That is doable, but requires significantly more computational effort and some experience with numerical methods to guarantee convergence and stability.

If you go that latter route, you might as well go with a proper nonlinear filter model, because that will need very similar methods to solve.