To further summarize, I want to create a function in matlab that finds the time domain signal $y(n)$ and its $n$ time components ($n=0,1,2,...$) given the numerator and denominator of a transfer function (filter) and the input sequence. I want to know if you can use the DFT and circular convolution to convolve $H(\omega)$ and $X(\omega)$ if the former is an IIR filter (there's a feedback).
You can't apply an IIR filter in the frequency domain, at least not without some approximation. As the name implies, the impulse response of an IIR filter is of infinite length which also means that you need infinitely high resolution in the frequency domain. Any DFT based implementation would require a finite frequency resolution, otherwise your FFT length becomes infinite.
Choosing a finite FFT length is the equivalent to truncating the impulse response.
In practice, of course, you can always come up with a good-enough approximation that's close enough to meet your specific requirements. Most IIRs have exponential decay so it dies off pretty fast. However, why would you? In terms of computational efficiency, memory consumption and latency, a direct IIR implementation will almost always beat any frequency domain algorithm by a wide margin. What's the problem you are trying to solve ?