Brief description of a FORTRAN 77 program is presented for calculating with the given accuracy eigenvalues, eigenfunctions and their first derivatives with respect to the parameter of the coupled parametric self-adjoined elliptic differential equations with the Dirichlet and/or Neumann type boundary conditions on the finite interval. The original problem is projected to the parametric homogeneous and nonhomogeneous 1D boundary-value problems for a set of ordinary second order differential equations which is solved by the finite element method. The program calculates also potential matrix elements – integrals of the eigenfunctions multiplied by their first derivatives with respect to the parameter. Parametric eigenvalues (so-called potential curves) and matrix elements computed by the POTHEA program can be used for solving the bound state and multi-channel scattering problems for a system of the coupled second-order ordinary differential equations with the help of the KANTBP programs. As a test desk, the program is applied to the calculation of the potential curves and matrix elements of Schr¨odinger equation for a system of three charged particles with zero total angular momentum.