Multiple solutions are common in various non-convex problems arising from industrial and scientific computing. Nonetheless, understanding the nontrivial solutions’ qualitative properties seems limited, partially due to the lack of efficient and reliable numerical methods. In this paper, we design a dedicated numerical method to explore these nontrivial solutions further. We first design an improved adaptive orthogonal basis deflation method by combining the adaptive orthogonal basis method with a bisection-deflation algorithm. We then apply the proposed new method to study the impact of domain changes on multiple solutions of certain nonlinear elliptic equations. When the domain varies from a circular disk to an elliptical disk, the corresponding functional value changes dramatically for some particular solutions, which indicates that these nontrivial solutions in the circular domain may become unstable in the elliptical domain. Moreover, several theoretical results on multiple solutions in the existing literature are verified. For the nonlinear sine-Gordon equation with parameter $λ,$ nontrivial solutions are found for $λ > λ_2,$ here $λ_2$ is the second eigenvalue of the corresponding linear eigenvalue problem. For the singularly perturbed Ginzburg-Landau equation, highly concentrated solutions are numerically verified, suggesting that their convergent limit is a delta function when the perturbation parameter goes to zero.