Construction solutions for Neumann problem with Hénon term in $ \mathbb{R}^2 $

In this paper, we are interested in the following boundary value problem

$ \begin{eqnarray*} \left\{ \begin{array}{ll} -\Delta u +u = \lambda |x-q_1|^{2\alpha_1}\cdots |x-q_n|^{2\alpha_n} u^{p-1}e^{u^p},\ \ u>0,\ \ \ & {\rm in}\ \Omega;\\ \frac{\partial u}{\partial\nu} = 0\ \ \ & {\rm on}\ \partial\Omega, \end{array} \right. \end{eqnarray*} $

where $ \Omega $ is a bounded domain in $ \mathbb{R}^2 $ with smooth boundary, points $ q_1,\ldots,q_n\in \Omega $, $ \alpha_1,\cdots,\alpha_n\in(0,\infty)\backslash\mathbb{N} $, $ \lambda>0 $ is a small parameter, $ 0< p <2 $, and $ \nu $ denotes the outer normal vector to $ \partial\Omega $. We construct solutions of this problem with $ k $ interior bubbling points and $ l $ boundary bubbling points, for any $ k\geq 1 $ and $ l\geq 1 $.