We study the following boundary value problem 0.1 $$\begin{aligned} \left\{ \begin{array}{ll} -div(a(x)\nabla u)+ a(x)u=0,\ u>0,\ \ \ \ &{}\quad \mathrm{in}\ \Omega ;\\ \frac{\partial u}{\partial \nu }=\lambda u^{p-1}e^{u^p},\ \ \ &{}\quad \mathrm{on}\ \partial \Omega , \end{array} \right. \end{aligned}$$ - d i v ( a ( x ) ∇ u ) + a ( x ) u = 0 , u > 0 , in Ω ; ∂ u ∂ ν = λ u p - 1 e u p , on ∂ Ω , where $$\Omega $$ Ω is a bounded domain in $$\mathbb {R}^2$$ R 2 with smooth boundary, $$\nu $$ ν is the unit outer normal vector of $$\partial \Omega $$ ∂ Ω , and a(x) is a positive smooth function in $${\overline{\Omega }}$$ Ω ¯ , $$\lambda >0$$ λ > 0 is a small parameter and $$0< p <2$$ 0 < p < 2 . We construct bubbling solutions to problem (0.1), which blow-up at points near a local maximum of a(x) on $$\partial \Omega $$ ∂ Ω .