In this paper, based on the bifurcation theory, we numerically calculate and visualize non-trivial multiple solutions of the singularly perturbed Neumann problem on a square using the Jacobi pseudospectral method. Starting from the non-zero solution of the corresponding problem, we numerically obtain multiple positive solutions with various symmetries by the branch switching method. The bifurcation diagrams are drawn, displaying the symmetry breaking bifurcation phenomena of the singularly perturbed Neumann problem. Numerical results demonstrate the effectiveness of this method.