Given a smooth bounded domain $\Omega$ in $\mathbb{ R}^2$, we study the following anisotropic Neumann problem $$ \begin{ cases} -\nabla(a(x)\nabla u)+a(x)u=\lambda a(x) u^{ p-1}e^{ u^p},\,\,\,\, u>0\,\,\,\,\, \textrm{ in}\,\,\,\,\, \Omega,\\[2mm] \frac{ \partial u}{ \partial\nu}=0\,\, \qquad\quad\qquad\qquad\qquad\qquad\qquad \ \ \ \ \,\qquad\quad\, \textrm{ on}\,\,\, \partial\Omega, \end{ cases} $$ where $\lambda>0$ is a small parameter, $0<p<2$, $a(x)$ is a positive smooth function over $\overline{ \Omega}$ and $\nu$ denotes the outer unit normal vector to $\partial\Omega$. Under suitable assumptions on anisotropic coefficient $a(x)$, we construct solutions of this problem with arbitrarily many mixed interior and boundary bubbles which concentrate at totally different strict local maximum or minimal boundary points of $a(x)$ restricted to $\partial\Omega$, or accumulate to the same strict local maximum boundary point of $a(x)$ over $\overline{ \Omega}$ as $\lambda\rightarrow0$.
Categories: math.AP