The effective and stable algorithms for numerical solution with the given accuracy of the parametrical two-dimensional (2D) boundary value problem (BVP)are presented. This BVP formulated for self-adjoined elliptic differential equations with the Dirichlet and/or Neumann type boundary conditions on a finite region of two variables. The original problem is reduced to the parametric homogeneous 1D BVP for a set of ordinary second order differential equations (ODEs). This reduction is implemented by using expansion of the required solution over an appropriate set of orthogonal eigenfunctions of an auxiliary Sturm-Liouville problem by one of the variables. Derivatives with respect to the parameter of eigenvalues and the corresponding vector-eigenfunctions of the reduced problem are determined as solutions of the parametric inhomogeneous 1D BVP. It is obtained by taking a derivative of the reduced problem with respect to the parameter. These problems are solved by the finite-element method with automatical shift of the spectrum. The presented algorithm implemented in Fortran 77 as the POTHEA program calculates with a given accuracy a set ∼ 50 of eigenvalues (potential curves), eigenfunctions and their first derivatives with respect to the parameter, and matrix elements that are integrals of the products of eigenfunctions and/or the derivatives of the eigenfunctions with respect to the parameter. The calculated potential curves and matrix elements can be used for forming the variable coefifcients matrixes of a system of ODEs which arises in the reduction of the 3D BVP (d = 3) in the framework of a coupled-channel adiabatic approach or the Kantorovich method. The efifciency and stability of the algorithm are demonstrated by numerical analysis of eigensolutions 2D BVP and evaluated matrix elements which apply to solve the 3D BVP for the Schro¨dinger equation in hyperspherical coordinates describing a Helium atom with zero angular momentum with help of KANTBP program.