Second-order conformal quantum superintegrable systems in 2 dimensions are Laplace equations on a manifold, with an added scalar potential and 3 independent 2nd order conformal symmetry operators. They encode all the information about 2D Helmholtz or time-independent Schroedinger superintegrable systems in an efficient manner: Each of these systems admits a quadratic symmetry algebra (not usually a Lie algebra) and is multiseparable. The separation equations comprise all of the various types of hypergeometric and Heun equations in full generality. In particular, they yield all of the 1D Schroedinger exactly solvable (ES) and quasi-exactly solvable (QES) systems related to the Heun operator. The separable solutions of these equations are the special functions of mathematical physics. The different systems are related by Staeckel transforms, by the symmetry algebras and by Bocher contractions of the conformal algebra so(4,C) to itself, which enables all systems to be derived from a single one: the generic potential on the complex 2-sphere. Distinct separable bases for a single Laplace system are related by interbasis expansion coefficients which are themselves special functions, such as the Wilson polynomials. Applying Bocher contractions to expansion coefficients for ES systems one can derive the Askey scheme for hypergeometric orthogonal polynomials. This approach facilitates a unified view of special function theory, incorporating hypergeometric and Heun functions in full generality.