Artificial compressibility methods aim to reduce the stiffness of the compressible Navier-Stokes equations by artificially decreasing the velocity of acoustic waves in the fluid. This approach has originally been developed as an alternative to the incompressible Navier-Stokes equations as this avoids the resolution of a Poisson equation. This paper extends the method to anisothermal low Mach number flows, allowing the simulations of subsonic flows submitted to large temperature variations, including dilatational effects. The procedure is shown to be stable and accurate using a finite-difference method in a staggered grid system for the simulation of strongly anisothermal turbulent channel flow. The highly scalable nature of the approach is well suited to complex high-fidelity simulations and GPU processing.