The photon statistics and, moreover, the density matrix (quantum state) of a single light mode can be sampled using homodyne detection. That is, the density matrix is computed by averaging a set of sampling functions with respect to the measured quadrature values. We develop a practical procedure for evaluating these functions. The algorithm is simple, stable, has small computer-memory requirements, and is fast enough for real-time data processing. (Interestingly, our method involves unnormalizable solutions of the stationary Schrodinger equation. We develop the annihilation-and-creation formalism for these solutions and derive their semiclassical approximations.) We quantify the bin width required to determine the density matrix up to a maximal quantum number. Finally, we show how the statistical errors of the reconstructed density matrix can be determined from the measured data.