In this paper, we first compute the multiple non-trivial solutions of the Schrödinger equation on a square, by using the Liapunov-Schmidt reduction and symmetry-breaking bifurcation theory, combined with Legendre pseudospectral methods. Then, starting from the non-trivial solution branches of the corresponding nonlinear problem, we further obtain the whole positive solution branch with $D_4$ symmetry of the Schrödinger equation numerically by pseudo-arclength continuation algorithm. Next, we propose the extended systems, which can detect the fold and symmetry-breaking bifurcation points on the branch of the positive solutions with $D_4$ symmetry. We also compute the multiple positive solutions with various symmetries of the Schrödinger equation by the branch switching method based on the Liapunov-Schmidt reduction. Finally, the bifurcation diagrams are constructed, showing the symmetry/peak breaking phenomena of the Schrödinger equation. Numerical results demonstrate the effectiveness of these approaches.