  Relativistic two-photon decay rates of the $2s_{1/2}$ and $2p_{1/2}$ states
towards the $1s_{1/2}$ ground state of hydrogenic atoms are calculated by using
numerically exact energies and wave functions obtained from the Dirac equation
with the Lagrange-mesh method. This approach is an approximate variational
method taking the form of equations on a grid because of the use of a Gauss
quadrature approximation. Highly accurate values are obtained by a simple
calculation involving different meshes for the initial, final and intermediate
wave functions and for the calculation of matrix elements. The accuracy of the
results with a Coulomb potential is improved by several orders of magnitude in
comparison with benchmark values of the literature. The general requirement of
gauge invariance is also successfully tested, down to rounding errors. The
method provides high accuracies for two-photon decay rates of a particle in
other potentials and is applied to a hydrogen atom embedded in a Debye plasma
simulated by a Yukawa potential.
