Summary: The lattice Boltzmann method (LBM) is extended to allow stable and consistent thermal evolution. Two multi-speed schemes are investigated; a three-speed, 34-state system that contains an artifact in the energy evolution equation, and a four-speed, 52-state system that removes this artifact. A dual-rate collision operator is used to achieve variable Prandtl number (Pr). The schemes are tested on the decay of a sinusoidal thermal perturbation that excites only the purely decaying eigenmode of this fully coupled system. Results over a range of transport coefficients, velocities, Pr and the allowable temperature range agree with theory for both schemes. However, under most flow conditions, the schemes eventually go unstable, a previously noted problem with multi-speed thermal models. We stabilize the scheme by identifying a temperature-dependent factor in the equilibrium function that leads directly to the removal of Galilean-invariance artifact, and relax the requirement of instantaneous accuracy of this factor. This results in a stable scheme but introduces artificial thermal diffusion strongly dependent on the bulk velocity. Empirically, we develop a scheme to push this error to higher-order. The resultant multi-speed schemes are also used to simulate Couette flow with a temperature gradient, and produce results that agree well with theory.