Solutions to the Laplace equation are required in many branches of groundwater hydrology. The most difficult class of problems involve generating solutions to the Laplace equation when the domain is an irregular shape and/or has a free surface. In this case, solutions can be generated by numerical Methods, e.g., finite difference, finite element, or boundary integer element methods. However, due to the computational requirements of these methods, they are not often incorporated into “practical” models of subcatchment hydrology. This paper, illustrates an alternative method to finite difference and finite elements for generating solutions to saturated and partiallysaturated porous media flow problems with mixed and irregular boundaries. The solution strategy is based on approximating the boundary conditions with orthonormal sequence expansions and developing appropriate recurrence relations for the coefficients of the expansion. In practice, the solution methodology involves truncation of the sequence to meet a desired accuracy. In past research, analytical series solutions were limited to problems with simple domains. However, we show that this method with modifications, can also be applied to complex domains. We also show that the method is computationally very efficient, and rapid solutions can be generated on basic PC computers. Several examples will be discussed.