We describe a variational method to solve the Holstein model for an electro
n coupled to dynamical, quantum phonons on an infinite lattice. The variati
onal space can be systematically expanded to achieve high accuracy with mod
est computational resources (12-digit accuracy for the one-dimensional pola
ron energy at intermediate coupling). We compute ground-state and low-lying
excited-state properties of the model at continuous values of the wave vec
tor k in essentially all parameter regimes. Our results for the polaron ene
rgy band, effective mass, and correlation functions compare favorably with
those of other numerical techniques, including the density-matrix renormali
zation-group technique, the global-local method, and the exact diagonalizat
ion technique. We find a phase transition for the first excited state betwe
en a bound and unbound system of a polaron and an additional phonon excitat
ion. The phase transition is also treated in strong-coupling perturbation t
heory.