  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$.
