An artificial neural network (ANN) method was developed for solving boundary value problems (BVPs) on an arbitrary irregular domain in such a manner that all Dirichlet and/or Neuman boundary conditions (BCs) are automatically satisfied. Exact satisfaction of BCs is not available with traditional numerical solution techniques such as the finite element method (FEM). The ANN is trained by reducing error in the given differential equation (DE) at certain points within the domain. Selection of these points is significantly simpler than the often difficult definition of meshes for the FEM. The approximate solution is continuous and differentiable, and can be evaluated at any location in the domain independent of the set of points used for training. The continuous solution eliminates interpolation required of discrete solutions produced by the FEM. Reducing error in the DE at a particular location in the domain does not necessarily imply improvement in the approximate solution there. A theorem was developed, proving that the solution will improve whenever error in the DE is reduced at all locations in the domain during training. The actual training of ANNs reasonably approximates the assumptions required by the proof. This dissertation offers a significant contribution to the field by developing a method for solving BVPs where all BCs are automatically satisfied. It had already been established in the literature that such automatic BC satisfaction is beneficial when solving problems on rectangular domains, but this dissertation presents the first method applying the technique to irregular domain shapes. This was accomplished by developing an innovative length factor. Length factors ensure BC satisfaction extrapolate the values at Dirichlet boundaries into the domain, providing a solid starting point for ANN training to begin. The resulting method has been successful at solving even nonlinear and non-homogenous BVPs to accuracy sufficient for typical engineering applications.