How does MATLAB handle IIR filters?
This isn't really a MATLAB-specific issue; I see a couple more general questions:

How do you implement a digital IIR filter?

You can apply any general digital filter by convolving its impulse response with the signal that you want to filter. That looks like:

$$
y[n] = \sum_{k=0}^{N-1} x[k] h[n-k]
$$

This works great for FIR filters, but you run into complications with IIR systems, because their impulse responses are infinitely long ($N \rightarrow \infty$ in the above sum). That makes the filter output difficult to calculate using the straightforward convolution sum above. To actually implement an IIR filter, we need a way to express the filter output in a way that is more computationally tractable.

The attribute of an IIR filter that makes its impulse response infinitely long is the system's recursive nature; there is feedback from the output back to the filter input. This means that the output of an IIR filter at any given instant of time can depend upon two separate information sources:


The input signal to the filter, past and present, and
Past values of the filter output signal.


This leads to the difference equation representation of the IIR system:

$$
a_0 y[n] = \sum_{k=1}^{M-1} a_k y[n-k] + \sum_{k=0}^{N-1} b_k x[n-k]
$$

This representation illustrates that we can implement an IIR filter by calculating a weighted sum of $M$ past outputs of the filter and $N$ past inputs of the filter (in practice, $N$ and $M$ could and often are the same; the order of the filter when defined as above is $max(M,N) - 1$). This is a closed-form expression that captures the behavior of the filter fully, suitable for automated implementation, and is the equation referenced in the MATLAB help excerpt that you showed in your question. Note that $a_0$ is almost exclusively assumed to be equal to $1$; I only included it because it is shown in the quoted help text, referenced as a(1).

So to specifically answer your question, as the text suggests, MATLAB implements IIR systems using the above equation. When you design an IIR filter with MATLAB or otherwise, you get two sets of coefficients ($a_k$ and $b_k$ in the difference equation) that define how to weight past filter inputs and outputs in the filter implementation. There are some other subtleties that are sometimes relevant; for instance, the text refers to using a "direct-form II transposed" realization of the filter. You then apply the filter by calculating the difference equation above for each desired output sample.

There are in fact a number of topologies that can be used to implement digital filters, each with their own tradeoffs. DF2T is often used because of its more efficient structure (it contains a minimum number of delay elements) and its improved robustness to roundoff error (since the feedforward coefficients $b_k$ are implemented first).

How does MATLAB design Butterworth filters?

The documentation for the butter function does not specify how it generates the discrete approximation to the analog Butterworth filter prototype. However, by peeking into butter.m, you find:

% step 5: Use Bilinear transformation to find discrete equivalent:
if ~analog,
    [a,b,c,d] = bilinear(a,b,c,d,fs);
end


So it seems that it uses the bilinear transform to map the Butterworth filter prototype to a digital filter realization.