The simulation of turbulence in the boundary region of a tokamak is crucial for understanding and optimizing the performance of fusion reactors. Solving a kinetic model for neutral particles deterministically using the method of characteristics invokes the solution of integral equations. Upon discretization, this requires the solution of dense, computationally expensive linear systems. This work employs low-rank linear algebra techniques to address this issue. In particular, we use hierarchical matrix approximations to significantly reduce the computational cost of assembling and solving these linear systems, leading to substantial savings in both time and memory. The hierarchical matrix method is implemented and tested within the GBS simulation code for boundary simulations, achieving over 90% reduction in computation time and memory, and enabling simulations with unprecedented spatial resolution for neutral particles.