Particle-in-cell (PIC) methods allow the study of plasma behavior by computing the trajectories of finite-size particles under the action of external and self-consistent electric and magnetic fields defined in a grid of points. In this work, the Finite Element Method (FEM) is used in order to obtain the self-consistent fields. An electrostatic PIC-FEM computational code for simulation of one-dimensional (1D) and two-dimensional (2D) plasmas was developed based on two available and independent codes: the first one a 1D PIC code that uses the Finite Difference Method and the other a FEM code developed at the Instituto de Estudos Avancados (IEAv). The Poisson equation is solved and periodic boundary conditions are used. The ion background that neutralizes the total plasma charge is kept fixed and uniformly distributed in the domain of study. The code is tested by studying the fluctuations of the plasma in thermal equilibrium. In thermal equilibrium a plasma sustain fluctuations of various collective modes of electrostatic oscillations, whose spectral distribution can be analytically obtained by using the fluctuation-dissipation theorem and the Kramers-Kronig relation. In both 1D and 2D cases, there are excellent agreement between the spectral distribution curves predicted theoretically and those obtained by simulation for finite size particles and long wavelengths.