The optical properties of photonic crystal slabs are often simulated with general three-dimensional methods such as finite-difference time-domain. Here we develop a multipole modal method, which is specialized to exploit two symmetries of the photonic crystal slab: the slab's vertically invariant nature allows the field to be expressed in Bloch modes, while the cylindrical inclusions allow the Bloch modes themselves to be expressed in the multipole basis. We find the multipole method to be fast and efficient in finding the Bloch modes, with convergence approaching the numerical accuracy possible. By matching the Bloch modes to plane waves at the top and bottom interfaces of the photonic crystal, the field scattered by the slab is calculated. Values of transmittance and reflectance accurate to 2-3 digits are easily and quickly achieved, whereas 5-6 digits are possible with greater numbers of modes and plane waves in field expansions. Higher accuracy is limited by Gibbs-related phenomena arising from the matching at the interface of necessarily discontinuous Bloch modes to necessarily continuous plane waves. We believe this limit may be present in all modal methods that use Bloch modes to expand the field within the photonic crystal.