  We present a finite element toolbox for the computation of Bogoliubov-de
Gennes modes used to assess the linear stability of stationary solutions of the
Gross-Pitaevskii (GP) equation. Applications concern one (single GP equation)
or two-component (a system of coupled GP equations) Bose-Einstein condensates
in one, two and three dimensions of space. An implementation using the free
software FreeFem++ is distributed with this paper. For the computation of the
GP stationary (complex or real) solutions we use a Newton algorithm coupled
with a continuation method exploring the parameter space (the chemical
potential or the interaction constant). Bogoliubov-de Gennes equations are then
solved using dedicated libraries for the associated eigenvalue problem. Mesh
adaptivity is proved to considerably reduce the computational time for cases
implying complex vortex states. Programs are validated through comparisons with
known theoretical results for simple cases and numerical results reported in
the literature.
