We propose new computational schemes and algorithms of the finite element method for solving elliptic multidimensional boundary value problems with variable coefficients at derivatives in a polyhedral d-dimensional domain, aimed at describing collective models of atomic nuclei. The desired solution is sought in the form of an expansion in the basis of piecewise polynomial functions constructed in an analytical form by joining Hermite interpolation polynomials and their derivatives on the boundaries of neighboring finite elements having the form of d-dimensional parallelepipeds. Calculations of the spectrum, quadrupole momentum and electric transitions of standard boundary value problems for the geometric collective model of atomic nuclei are analyzed.