Context. The study of the continuum radiative transfer problem inside circumstellar envelopes is both a theoretical and numerical challenge, especially in the frequency-dependent and multi-dimensional case. While approximate methods are easier to handle numerically, they often fail to accurately describe the radiation field inside complex geometries. For these cases, it is necessary to directly solve the radiative transfer equation numerically.Aims. We investigate the accuracy of the discontinuous Galerkin finite element method (DGFEM hereafter) applied to the frequency-dependent two-dimensional radiative transfer problem, and coupled with the radiative equilibrium equation. We next used this method in the context of axis-symmetric circumstellar envelopes.Methods. The DGFEM is a variant of finite element methods. It employs discontinuous elements and flux integrals along their boundaries, ensuring local flux conservation. However, as opposed to the classical finite element methods, the solution is discontinuous across element edges. We implemented this approach in a code and tested its accuracy by comparing our results with the benchmarks from the literature.Results. For all the tested cases, the temperatures profiles agree within one percent. Additionally, the emerging spectral energy distributions (SEDs) and images, obtained by ray-tracing techniques from the DGFEM emissivity, agree on average within 5% and 10%, respectively.Conclusions. We show that the DGFEM can accurately describe the continuum radiative transfer problem inside axis-symmetric circumstellar envelopes. Consecutively the emerging SEDs and images are also well reproduced. The DGFEM provides an alternative method (other than Monte-Carlo methods for instance) for solving the radiative transfer equation, and it could be used in cases that are more difficult to handle with the other methods.