A method is proposed for analysis of the gamma-ray pulse-height distribution (PHD) obtained with a Ge(Li) detector, which makes use of a fitting technique based on a nonlinear least-squares method. The successive steps of analysis are thus made from peak search to estimations of gamma-ray peak areas, peak locations, and their errors. Gamma-ray peaks and background continuum are expressed in a Gaussian distribution function; an exponential function is connected smoothly with the Gaussian at a certain point on the low-energy side of a gamma-ray peak, a newly defined Valley Count Function (VCF) expressing the flat-background continuum located between Compton edge and centroid of peak, and a quadratic polynomial expressing the ordinary background continuum. From the results of the analysis obtained with a computational code, it was found that, to obtain meaningful solutions in gamma-ray PHD analysis, the full width at half maximum (fwhm) must contain three or more channels, and the separation between the peaks should be larger than half the fwhm in the analysis of a composed PHD having two peaks. The validity of the method is shown by an analysis of the gamma-ray spectrum of Ta isotopes.