In plasma physics, the Boltzmann equation, which describes the evolution of the plasma over time, has a nonlinear term representing the collisions of various species of the plasma. Current plasma edge simulations do not take this collision effect into account, because of the difficulties in the accurate evaluation of this term. Using the Rosenbluth potential formalism, the collision operator can be written in terms of solutions of a Poisson and a biharmonic free space PDE. Due to the inherent axisymmetry of the input data, cylindrical coordinate solvers are preferred for efficient computation. Standard numerical techniques (based typically on finite differences and finite element approximations) encounter difficulties in achieving high order accuracy, especially in the computation of derivatives of the solution (required in the collision operator formulation), and in imposing radiation conditions at infinity. Our new solver achieves arbitrary order accuracy in cylindrical coordinates based on a combination of separation of variables, Fourier analysis and the careful solution of the resulting radial ODE. A weak singularity arises in the the continuous Fourier transform of the solution that can be handled effectively with special purpose quadrature rules and spectral accuracy can be achieved in derivatives without loss of precision. This is joint work with Leslie Greengard.