Among the classical solution methods for partial differential equations, the method of separation of variables assumes that, on elementary domains such as rectangles (rectangular boxes) or disks (balls), the solution can be represented as a product of functions, each depending only on a single coordinate variable. One then determines the separated functions and constructs the general solution by the principle of superposition. In this work, we address whether such an approach is possible in more general domains. Instead of separating the solution with respect to variables, we propose separating the differential operator itself into coordinate-direction differential operators. With this approach, for general elliptic partial differential equations with variable coefficients, the operator is decomposed into differential operators corresponding to each coordinate direction. By employing one-dimensional Green’s functions, one obtains a system of one-dimensional integral equations whose number equals the spatial dimension. During this sequence of formulations, new auxiliary functions arise and must be determined. In this talk, we investigate the theoretical properties of numerical solutions obtained by discretizing the resulting integral equations. In particular, we focus on their connection with finite difference methods, various theoretical properties, stability of the solutions, and error analysis.