Computational schemes of the Galerkin type method (GTM) and finite elements method (FEM) for solving elliptic multidimensional boundary value problems (BVPs) with variable coefficients of derivatives in a polyhedral d-dimensional domain, aimed at describing collective models of atomic nuclei are presented. The solution is sought in the form of an expansion in the GTM basis and/or in the FEM basis of piecewise polynomial functions constructed in analytical form by joining Hermite interpolation polynomials and their derivatives at the boundaries of neighboring finite elements, which have the form of d-dimensional parallelepipeds. The BVPs are formulated and analyzed for collective models including the mixed derivative of the two-dimensional vibrational part of the five-dimensional Hamiltonian in the representation of the nuclear spin angular momentum in the intrinsic reference frame defined by three Euler angles. Benchmark calculations demonstrate performance and robustness of the approach when applied to calculate the lower part of the energy spectrum and the reduced electric transition probabilities in quadrupole collective models of atomic nuclei. The calculations of the band spectrum of Gd isotope using tabulated variable coefficients of the BVP evaluated in the self-consistent relativistic mean-field model revealed a possibility of quasicrossing of energy levels belonging to different rotational bands of a nucleus at high spin values. © The Author(s), under exclusive license to Springer Nature Switzerland AG 2024.