Reliable global elucidation of (subsets of) self-consistent field solutions is required for continued development and application of computational approaches that utilize these solutions as reference wavefunctions. We report the derivation and implementation of a stochastic approach to perform global elucidation of self-consistent field solutions by exploiting the connection between global optimization and global elucidation problems. We discuss the design of the algorithm through combining basin-hopping search algorithms with a Lie algebraic approach to linearize self-consistent field solution space, while also allowing preservation of desired spin-symmetry properties of the wavefunction. The performance of the algorithm is demonstrated on minimal basis C2v H4 due to its use as a model system for global self-consistent field solution exploration algorithms. Subsequently, we show that the model is capable of successfully identifying low-lying self-consistent solutions of benzene and NO2 with polarized double-zeta and triple-zeta basis sets and examine the properties of these solutions.