The design optimization and analysis of charged particle beam systems employing intense beams requires a robust and accurate Poisson solver. This paper presents a new type of Poisson solver which allows the effects of space charge to be elegantly included into the system dynamics. This is done by casting the charge distribution function into a series of basis functions, which are then integrated with an appropriate Green's function to find a Taylor series of the potential at a given point within the desired distribution region. In order to avoid singularities, a Duffy transformation is applied, which allows singularity-free integration and maximized convergence region when performed with the help of Differential Algebraic methods. The method is shown to perform well on the examples studied. Practical implementation choices and some of their limitations are also explored.