We present a coherent and effective theoretical framework to systematically construct numerically exact nonlinear solitary waves from their respective linear limits. First, all possible linear degenerate sets are classified for a harmonic potential using lattice planes. For a generic linear degenerate set, distinct wave patterns are identified in the near-linear regime using a random searching algorithm by suitably mixing the linear degenerate states, followed by a numerical continuation in the chemical potential extending the waves into the Thomas–Fermi regime. The method is applied to the two-dimensional, one-component Bose–Einstein condensates, yielding a spectacular set of waveforms. Our method opens a remarkably large program, and many more solitary waves are expected. Finally, the method can be readily generalized to three dimensions, and also multi-component condensates, providing a highly powerful technique for investigating solitary waves in future works.