The one-dimensional Helmholtz equation, ε2 u xx −u=f(x), arises in many applications, often as a component of three-dimensional fluids codes. Unfortunately, it is difficult to solve for ε≪1 because the homogeneous solutions are exp (±x/ε), which have boundary layers of thickness O(1/ε). By analyzing the asymptotic Chebyshev coefficients of exponentials, we rederive the Orszag–Israeli rule [16] that $$N \approx {3 \mathord{\left/ {\vphantom {3 {\sqrt \varepsilon }}} \right. \kern-\nulldelimiterspace} {\sqrt \varepsilon }}$$ Chebyshev polynomials are needed to obtain an accuracy of 1% or better for the homogeneous solutions. (Interestingly, this is identical with the boundary layer rule-of-thumb in [5], which was derived for singular functions like tanh([x−1]/ε).) Two strategies for small ε are described. The first is the method of multiple scales, which is very general, and applies to variable coefficient differential equations, too. The second, when f(x) is a polynomial, is to compute an exact particular integral of the Helmholtz equation as a polynomial of the same degree in the form of a Chebyshev series by solving triangular pentadiagonal systems. This can be combined with the analytic homogeneous solutions to synthesize the general solution. However, the multiple scales method is more efficient than the Chebyshev algorithm when ε is very, very tiny.