An algorithm, which reduces to velocity Verlet in the limit of zero friction, is obtained for the generalized Langevin equation. The formulation presented is unique in that the velocities are based on a direct second order Taylor expansion of the inertial forces. The resulting finite difference equation is compared with a previous formulation of Verlet-based Langevin dynamics. The equations are implemented to study the physical properties of dense neon and liquid water at constant temperatures as a function of the friction rate γ. The results show canonical ensembles can be obtained at moderate friction rates without recurring to switching functions to control energy losses, or overdamping. The LD methodology is further applied to a protein in solution and compared to MD simulations which use a more conventional scaling function for temperature control. Analyses of the MD trajectories show the protein remains at either lower or higher temperature compared to the solvent, depending upon the treatment of the potential function at the system's boundary. On the contrary, LD simulations show a rapid and even equilibration between protein and solvent in all cases.