This study presents a semi-analytical approach to solve the Laplace equation in arbitrarily shaped two-dimensional domains. The method is meshless and addresses boundary value problems that include both pure Dirichlet and mixed Dirichlet-Neumann boundary conditions. The solution is obtained via a weighted superposition of harmonic polynomials which are obtained via an orthonormalization approach. We show that the approach is efficient in terms of number of operations. The numerical solution is convergent and exact (within machine precision) given a sufficient number of terms in the series. Moreover, the method offers several advantages over traditional approaches. Advantages include providing analytical expressions for the stream function and velocity components when solving potential flow problems. There are also important implications for the storage of model results; the method offers extremely low cost data storage. In the paper, several example applications involving arbitrary domains are presented. The results obtained are compared with known analytical solutions.