We consider the calculation schemes for solving elliptic boundary-value problems (BVPs) within the framework of the Kantorovich method that provides the reduction of an elliptic BVP to a system of coupled second-order ordinary differential equations (ODEs). The surface basis functions of the expansion depend on the independent variable of the ODEs parametrically. Here we use the basis functions calculated by means of the finite element method(FEM), as well as the probe parametric surface basis functions calculated in the analytical form. We propose new calculation schemes and algorithms for solving the parametric self-adjoint elliptic boundary-value problem (BVP) in a 2D finite domain, using high-accuracy finite element method (FEM) with rectangular and triangular elements. The algorithm and the programs calculate with the given accuracy the eigenvalues, the surface eigenfunctions and their first derivatives with respect to the parameter of the BVP for parametric self-adjoint elliptic differential equation with the Dirichlet and/or Neumann type boundary conditions on the 2D finite domain, and the potential matrix elements, expressed as integrals of the products of surface eigenfunctions and/or their first derivatives with respect to the parameter. The parametric eigenvalues (potential curves) and the potential matrix elements computed by the program can be used for solving bound-state and multi-channel scattering problems for systems of coupled second-order ODEs by means of the Kantorovich method. We demonstrate the efficiency of the proposed calculation schemes and algorithms in benchmark calculations of 2D elliptic BVPs describing quadrupole vibrations of a collective nuclear model.