We propose new symbolic-numerical algorithms implemented in Maple-Fortran environment for solving the parametric self-adjoint elliptic boundary-value problem (BVP) in a 2D finite domain, using high-accuracy finite element method (FEM) with triangular elements and high-order fully symmetric Gaussian quadratures with positive weights, and no points are outside the triangle (PI type). The algorithms 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. We demonstrated an efficiency of algorithms and program by benchmark calculations of helium atom ground state. © 2017, Springer International Publishing AG.