  We describe a systematic approach for the efficient numerical solution of
nonlinear Schr\"odinger-type partial differential equations of the form $(K +V
+ g|\psi|^2)\psi=0$, with an energy operator $K$, a scalar potential $V$, and a
scalar parameter $g$. Instrumental to the approach are developments in
numerical linear and nonlinear algebra, specifically numerical parameter
continuation. We demonstrate how a continuous sequence of solutions can be
obtained regardless of their stability, so that finally the spectrum of stable
and unstable solutions in the specified parameter range is fully revealed. The
method is demonstrated for the GL equation in a three-dimensional
superconducting domain with an inhomogeneous magnetic field, a numerically
demanding problem known to have an involved solution landscape.
