A numerically efficient Fredholm formulation of the billiard problem is presented. The standard solution in the framework of the boundary integral method in terms of a search for roots of a secular determinant is reviewed first. We next reformulate the singularity condition in terms of a flow in the space of an auxiliary one-parameter family of eigenproblems and argue that the eigenvalues and eigenfunctions are analytic functions within a certain domain. Based on this analytic behavior, we present a numerical algorithm to compute a range of billiard eigenvalues and associated eigenvectors by only two diagonalizations.