Really? Can you expand on what parts are tricky to get right? I would have thought for simple 2D problems like this (with simple geometry and where O(N^3) methods are fine so that you don’t need iterative solvers or FMM) you can just mesh the domain yourself, setup the integral equation and discretize it brutally. At least that seems like a fun learning exercise (although less flexible for more complicated geometries, I agree).
(Also, since the Green function for the Poisson equation in the disk is known analytically, you can possibly even not have to mesh the circle?)