Third- and fourth-order FEM schemes with multivariate Hermite interpolation polynomials of a d-dimensional hypercube for solving boundary value problems (BVPs) on hyperparallelepipedal meshes are elaborated. An exactly solvable model of a system of several identical particles with pair oscillator interaction known as the Moshinsky atom is used as a test example. To describe the energy spectra of symmetric and antisymmetric bound states, the 2-,3-,4-, and 5- dimensional BVPs with Dirichlet and Neumann boundary conditions on a nonrectangular domain are formulated. To generate new FEM schemes with mixed partial derivatives, additional affine coordinate transformations are applied. Benchmark calculations of the BVPs confirm the order of declared FEM schemes.