Let \(\varepsilon >0\) small enough and \(\lambda>0\) large. Assume that \(\Omega\subset {\mathbb R}^2\) is a bounded domain with smooth boundary. The authors construct solutions to the Liouville equation \[ \Delta u+\varepsilon^2e^u=0\qquad \text{in}\;\Omega\,, \] under the Robin boundary condition \[ \frac{\partial u}{\partial\nu}+\lambda u=0\qquad \text{on}\;\partial\Omega\,. \] The solutions constructed exhibit concentration as \(\varepsilon\rightarrow 0\) and simultaneously as \(\lambda\rightarrow +\infty\) at points that get close to the boundary. The authors show that, in general, the set of solutions of this problems exhibits a richer structure than the problem with Dirichlet boundary conditions. The proof combines arguments for a projected version of the nonlinear equation with a careful expansion of the energy of the ansatz and estimates of the Green function.