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 from 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.