IIR filters are also linear (and time-invariant). what that means is that the mapping from input samples $x[n]$ to output samples $y[n]$ satisfies the superposition (and time-shifting) properties:
if $$y[n] \triangleq \mbox{LTI} \{ x[n] \}$$
then $$y_1[n] + y_2[n] = \mbox{LTI} \{ x_1[n] + x_2[n] \}$$
and $$ y[n-D] = \mbox{LTI} \{ x[n-D] \} $$
and, as a result, the convolution summation:
$$y[n] = \mbox{LTI} \{ x[n] \} = \sum_{k=-\infty}^{+\infty} h[k] x[n-k]$$
for some impulse response $h[n]$ (which are the coefficients for the FIR) that is dependent only on the filter, not the input. that's what linear and time-invariant means.
but in both cases, FIR or IIR, the mapping of the specifications of the filter to those coefficients $h[n]$ is the "filter design" process. why would that be linear? there is no reason for that.
so, perhaps one method might be you make a guess at what the coefficients are, see what you get regarding the coefficients, and then somehow tweek or adjust your guess to be a better guess in the next recursion. now algorithms like the Parks-McClellan (based on the Remez exchange algorithm) are more complicated than "guess and tweek", but that is what they are doing.
you can get very good results with the "windowing method" (if you use a good window function, like the Kaiser) and that is not recursive.