  We propose a First-Order System Least Squares (FOSLS) method based on
deep-learning for numerically solving second-order elliptic PDEs. The method we
propose is capable of dealing with either variational and non-variational
problems, and because of its meshless nature, it can also deal with problems
posed in high-dimensional domains. We prove the $\Gamma$-convergence of the
neural network approximation towards the solution of the continuous problem,
and extend the convergence proof to some well-known related methods. Finally,
we present several numerical examples illustrating the performance of our
discretization.
