Electrostatic structures have been observed in many regions of space plasmas, including the solar wind, the magnetosphere, the auroral acceleration region, and in association with shocks, turbulence, and magnetic reconnection. Due to potentially large amplitude of electric fields within these structures, their effects on particle heating, scattering, or acceleration can be important. One possible theoretical description of some of these structures is the concept of Bernstein-Greene-Kruskal (BGK) modes, which are exact nonlinear solutions of the Vlasov-Poisson system of equations in collisionless kinetic theory. BGK modes have been studied extensively for many decades, predominately in one dimension (1D), although there have been observations showing that some of these structures have clear 3D features. While there have been approximate solutions of higher dimensional BGK modes, an exact 3D BGK mode solution in a finite magnetic field has not been found yet. Recently we have constructed exact solutions of 2D BGK modes in a magnetized plasma with finite magnetic field strength in order to gain insights of the ultimate 3D theory [Ng, Bhattacharjee, and Skiff, Phys. Plasmas 13, 055903 (2006)]. Based on the analytic form of these solutions, as well as Particle-in-Cell (PIC) simulations, we will present numerical studies of their stability for different levels of background magnetic field strength.