A finite element approach to solve numerically the Takagi-Taupin equations expressed in a weak form is presented and applied to simulate X-ray reflectivity curves, spatial intensity distributions and focusing properties of bent perfect crystals in symmetric reflection geometry. The proposed framework encompasses a new formulation of the Takagi-Taupin equations, which appears to be promising in terms of robustness and stability and supports the Fresnel propagation of the diffracted waves. The presented method is very flexible and has the potential of dealing with dynamical X-ray or neutron diffraction problems related to crystals of arbitrary shapes and deformations. The reference implementation based on the commercial COMSOL Multiphysics software package is available to the relevant user community.