When two or more branches of a function merge, the Chebyshev series of u(λ) will converge very poorly with coefficients a n of T n (λ) falling as O(1/n α ) for some small positive exponent α. However, as shown in [J.P. Boyd, Chebyshev polynomial expansions for simultaneous approximation of two branches of a function with application to the one-dimensional Bratu equation, Appl. Math. Comput. 143 (2002) 189-200], it is possible to obtain approximations that converge exponentially fast in n. If the roots that merge are denoted as u 1 (λ) and u 2 (λ), then both branches can be written without approximation as the roots of (u − u 1 (λ))(u − u 2 (λ)) = u 2 + β(λ)u + γ (λ). By expanding the nonsingular coefficients of the quadratic, β(λ) and γ (λ), as Chebyshev series and then applying the usual roots-of-a-quadratic formula, we can approximate both branches simultaneously with error that decreases proportional to exp(−σ N ) for some constant σ > 0 where N is the truncation of the Chebyshev series. This is dubbed the "Chebyshev-Shafer" or "Chebyshev-Hermite-Padé" method because it substitutes Chebyshev series for power series in the generalized Padé approximants known variously as "Shafer" or "Hermite-Padé" approximants. Here we extend these ideas. First, we explore square roots with branches that are both real-valued and complex-valued in the domain of interest, illustrated by meteorological baroclinic instability. Second, we illustrate triply branched functions via roots of the Kepler equation, f (u; λ, ) ≡ u − sin(u) − λ = 0. Only one of the merging roots is real-valued and the root depends on two parameters (λ, ) rather than one. Nonetheless, the Chebyshev-Hermite-Padé scheme is successful over the whole two-dimensional parameter plane. We also discuss how to cope with poles and logarithmic singularities that arise in our examples at the extremes of the expansion domain.