Presently, state‐of‐the‐art theoretical studies of elementary excitations at surfaces involve large, time consuming, and often cumbersome calculations, which are not very practical for surfaces with low symmetry. A new method has just been developed that directly yields more information and is computationally faster (by an order of magnitude) than any current method available. This method will allow calculations of complicated surface systems not previously possible. The ingredients of the theory are any semi‐infinite system of atomic layers, any set of basis functions spanning each layer, and one single matrix which depends only on the interaction matrix elements among the atomic layers. What is remarkable is that the eigenanalysis of this matrix gives a complete quantum mechanical description of the elementary excitations of the system. Therefore with just a diagonalization of this matrix, the knowledge of its eigenvalues and eigenvectors give directly (1) the projected bandstructure (without having to diagonalize the bulk Hamiltonian for a large number of K⊥’s); (2) the bona fide surface bands with high accuracy (without the need of a bulk Green function first); (3) the surface Green function as a function of energy E+iδ, in the correct limit that δ→o! (without the need to solve slowly convergent iterative equations for complexenergies); and (4) the precise decay length, orbital character, and symmetry of the wave function. The method, although versatile, is also conceptually simple. All the properties mentioned above are obtained from a simple set of rules. A tutorial discussion of the formalism will be presented along with several instructive examples illustrating its use.