We describe a variational method to solve the Holstein model for an electron coupled to dynamical, quantum phonons on an infinite lattice. The variational space can be systematically expanded to achieve high accuracy with modest computational resources (12-digit accuracy for the 1d polaron energy at intermediate coupling). We compute ground and low-lying excited state properties of the model at continuous values of the wavevector $k$ in essentially all parameter regimes. Our results for the polaron energy band, effective mass and correlation functions compare favorably with those of other numerical techniques including DMRG, Global Local and exact diagonalization. We find a phase transition for the first excited state between a bound and unbound system of a polaron and an additional phonon excitation. The phase transition is also treated in strong coupling perturbation theory.