In this paper, we are concerned with stationary solutions to the following Gierer–Meinhardt system with saturation and source term under the homogeneous Neumann boundary condition: { A t = ε 2 Δ A − A + A 2 H ( 1 + k A 2 ) + σ 0 in Ω × ( 0 , ∞ ) , τ H t = D Δ H − H + A 2 in Ω × ( 0 , ∞ ) . Here, ε > 0 , τ ≥ 0 , k ≥ 0 , and Ω is a bounded smooth domain in R N . In this paper, we suppose Ω is an x N -axially symmetric domain and σ 0 is an x N -axially symmetric nonnegative function of class C α ( Ω ¯ ) , α ∈ ( 0 , 1 ) . For sufficiently small ε and sufficiently large D , we construct a multi-peak stationary solution peaked at arbitrarily chosen intersections of the x N -axis and ∂ Ω under the condition that 4 k ε − 2 N | Ω | 2 converges to some k 0 ∈ [ 0 , ∞ ) as ε → 0 . This extends the results of Kurata and Morimoto [K. Kurata, K. Morimoto, Construction and asymptotic behavior of the multi-peak solutions to the Gierer–Meinhardt system with saturation, Commun. Pure Appl. Anal. 7 (2008) 1443–1482] to the case σ 0 ( x ) ≥ 0 . Moreover, we study an effect of the source term σ 0 on a precise asymptotic behavior of the solution as ε → 0 .