The simplest example of a boundary value problem is the Dirichlet Poisson problem: we seek to recover a function, defined on a smooth domain, given its values on the boundary of the domain and the divergence of its gradient inside the domain. This problem has been studied for more than 200 years and has many applications in science and engineering. Analytic solutions are available only for a limited number of cases. Therefore, one has to use a numerical method. The goal in designing a numerical method is the ability to guarantee a number of correct digits in the numerical solution, in reasonable time, and in a black-box fashion. Surprisingly, such a robust, algorithmically scalable method does not exist for the Poisson problem. The Finite Element Method is the most popular tool for boundary value problems. FEM based approaches however, cannot operate in a black-box fashion. The main difficulties are unstructured mesh generation and ill-conditioning of the algebraic system of equations. In this talk I will review different approaches in solving the Poisson problem and present a new method based on a classical Fredholm integral equation formulation. The main components of the new method are a kernel-independent fast summation method, manifold surface representations, and superalgebraically accurate quadrature methods. This is join work with Lexing Ying of Caltech, and Denis Zorin of New York University.