First-order system least squares (FOSLS) is a methodology that offers an alternative to standard methods for solving partial differential equations. This paper studies the first-order system least-squares approach for scalar second-order elliptic boundary value problems with discontinuous coefficients. In a companion paper, ellipticity of an appropriately scaled least-squares bilinear form is established in a natural norm. For some geometries this ellipticity is independent of the size of the jumps in the coefficients. The occurrence of singularities at interface corners, cross-points, reentrant corners, and irregular boundary points is discussed, and a basis of singular functions with local support around singular points is established. This paper describes a method for including discrete versions of the singular basis functions together with standard finite element spaces in a least-squares format at little additional computational cost. The singular basis functions are constructed to match the jump conditions that arise at interfaces between regions of continuity of the diffusion coefficient. Because these basis functions must be approximated in practice, the resulting discretization is by nature nonconforming. This necessitates the establishment here of a general error estimate for FOSLS L2 minimization problems discretized by nonconforming finite elements. An advantage of the FOSLS formulation is that this estimate does not involve the consistency error term usually present in bounds for other methods. Based on this general estimate error bounds are derived for the finite element space that includes singular basis functions. Numerical tests are included that confirm these discretization error bounds. Finally, a multilevel method is developed for solving the discrete system that uses singular basis functions on all levels and its efficiency is demonstrated by the numerical results.