We implement in Maple and Mathematica an algorithm for constructing multivariate Hermitian interpolation polynomials (HIPs) inside a d-dimensional hypercube as a product of d pieces of one-dimensional HIPs of degree $$p'$$ in each variable, that are calculated analytically using the authors’ recurrence relations. The piecewise polynomial functions constructed from the HIPs have continuous derivatives and are used in implementations of the high-accuracy finite element method. The efficiency of our finite element schemes, algorithms and GCMFEM program implemented in Maple and Mathematica are demonstrated by solving reference boundary value problems (BVPs) for multidimensional harmonic and anharmonic oscillators used in the Geometric Collective Model (GCM) of atomic nuclei. The BVP for the GCM is reduced to the BVP for a system of ordinary differential equations, which is solved by the KANTBP 5 M program implemented in Maple.