In this paper, we present a novel approach utilizing the lattice Boltzmann method (LBM) to simulate non-Newtonian fluids characterized by the White-Metzner constitutive equation under low Reynolds number conditions. By recognizing that the elastic stress tensor has a non-vanishing trace, we are able to decompose it into its deviatoric and spherical components. We develop a theoretical reformulation of the White-Metzner model, where these components are computed separately. Specifically, our approach employs three distinct distribution functions: two scalar distributions for the evolution of momentum and elastic pressure, and one tensor distribution for the evolution of deviatoric stress. We provide a comprehensive derivation of the governing equations for both deviatoric stress and elastic pressure. Given that the elastic pressure provides insights into the elongation of polymers dissolved in the solvent, a four-roll mill geometry is utilized as a benchmark test to investigate the elastic instabilities that arise with increasing Weissenberg numbers. The flow behavior is examined numerically using a coupled lattice Boltzmann method, demonstrating the effectiveness of our approach in analyzing these instabilities.