We describe a special Gauss-Newton method for tracing solution manifolds with singularities of multiparameter systems. First we choose one of the parameters as the continuation parameter., and fix the others. Then we trace one-dimensional solution curves by using continuation methods. Singularities such as folds, simple and multiple bifurcations on each solution curve can be easily detected. Next, we choose an interval for the second continuation parameter, and trace one-dimensional solution curves for certain values in this interval. This constitutes a two-dimensional solution surface. The procedure can be generalized to trace a k-dimensional solution manifold. Numerical results in 1D, 2D and 3D second-order semilinear elliptic eigenvalue problems, given by Lions [1982] are reported.