Axial Green's function method is proposed to solve multi-dimensional second-order elliptic partial differential equations with the geometrical complexity of solution domains. Instead of directly dealing with 2D or 3D Green's functions, the proposed method systematically decomposes a multi-dimensional elliptic differential operator into 1D ones. Consequently, the method only requires 1D Green's function associated with each coordinate axis. Theoretical formulation and simple axial discretization result in the system of equations of which the solutions represent the numerical solution of the original multi-dimensional elliptic problem, which is otherwise intractable to obtain by the traditional multi-dimensional Green's function technique. The method employs analytical form of 1D Green's function and integral schemes so that it can potentially be more robust than both the finite difference method and the boundary element method. In addition, the simple discretization procedure, which uses only axially orthogonal straight lines and boundary data, demonstrates its convenience over the finite element method of comparable accuracy. Numerical examples verify these advantages by demonstrating the second-order convergence in various discretized models of complex geometries. AXIAL GREEN'S FUNCTION METHOD where the bracket implies [•] t= = lim t→ + •−lim t→ − •. (ii) Continuity: G(t, ) is continuous in I with respect to the variable t.