Relativistic dipolar to hexadecapolar polarizabilities of the ground state and some excited states of hydrogenic atoms are calculated by using numerically exact energies and wave functions obtained from the Dirac equation with the Lagrange-mesh method. This approach is an approximate variational method taking the form of equations on a grid because of the use of a Gauss quadrature approximation. The partial polarizabilities conserving the absolute value of the quantum number $kappa$ are also numerically exact with small numbers of mesh points. The ones where $|kappa|$ changes are very accurate when using three different meshes for the initial and final wave functions and for the calculation of matrix elements. The polarizabilities of the $n=2$ excited states of hydrogenic atoms are also studied with a separate treatment of the final states that are degenerate at the nonrelativistic approximation. The method provides high accuracies for polarizabilities of a particle in a Yukawa potential and is applied to a hydrogen atom embedded in a Debye plasma.