Based on the recently developed finite integration method for solving partial differential equations,we extend in this paper the advantages of the method in tackling singular perturbation problems with multiple boundary layers.Numerical results indicate that,even with the most simple numerical trapezoidal integration rule,the method provides very stable,efficient,and highly accurate solutions to the singular perturbation problems.Convergence error estimate of the method is derived.An adaptive refinement of integration points is also proposed for further improvement of the method in locating multiple boundary layers.Illustrative examples are given to compare with other existing numerical methods.