A new methodology is introduced for the numerical solution of Partial Differential Equations in general spatial domains. The methodology is based on the use of the well-known Alternating Direction Implicit (ADI) approach of Peaceman and Rachford in conjunction with one-dimensional and high-order accurate Fourier representations of non-periodic data, obtained by way of a certain "continuation method" introduced recently for the resolution of the Gibbs phenomenon. We construct a number of high-order convergent PDE solvers on the basis of this strategy. Unlike previous alternating direction methods for general domains of order higher than one, the new algorithms possess the desirable property of unconditional stability for general spatial domains; the computational time required for these methods to advance one time-step, in turn, grows in an essentially linear manner with the number of spatial discretization points. In particular, the new methodology yields significant advantages over traditional low-order methods for computations involving wave propagation and "large domains," as well as PDEs including diffusive terms. In all, we treat Dirichlet problems for the Heat Equation, the Poisson Equation, and the Wave Equation in two- and three-dimensional spatial domains with smooth boundaries. A stability analysis we present hinges upon the numerical evaluation of certain singular value decompositions. Numerical results arising from the implementations of two- and three-dimensional versions of the method for linear parabolic, hyperbolic and elliptic PDEs, exhibit unconditional stability and high-order convergence, in agreement with our theoretical results.