In this paper, we employ the weak Galerkin (WG) finite element method and the imaginary time method to compute both the ground state and the excited states in Bose-Einstein condensate (BEC) which is governed by the Gross-Pitaevskii equation (GPE). First, we use the imaginary time method for GPE to get the nonlinear parabolic partial differential equation. Subsequently, we apply the WG method to spatially discretize the parabolic equation. This yields a semi-discrete scheme, in which an energy function is explicitly defined. For the case β⩾0, we demonstrate that the energy is diminishing with respect to time t at each time step. Applying the backward Euler scheme for temporal discretization yields a fully discrete scheme. For the case β=0, we provide a mathematical justification, establishing the convergence analysis for the numerical solution of the ground state. Moreover, based on the theory of solving eigenvalue problems using the WG method, we present the error estimates between the ground state and its numerical solution under the H1 and L2 norms. Numerical experiments are provided to illustrate the effectiveness of the proposed schemes. Moreover, the results indicate that our method also can compute the first excited state, achieving optimal convergence orders.