  We develop a spectrally accurate numerical method to compute solutions of a
model partial differential equation used in plasma physics to describe
diffusion in velocity space due to Fokker-Planck collisions. The solution is
represented as a discrete and continuous superposition of normalizable and
non-normalizable eigenfunctions via the spectral transform associated with a
singular Sturm-Liouville operator. We present a new algorithm for computing the
spectral density function of the operator that uses Chebyshev polynomials to
extrapolate the value of the Titchmarsh-Weyl $m$-function from the complex
upper half-plane to the real axis. The eigenfunctions and density function are
rescaled and a new formula for the limiting value of the $m$-function is
derived to avoid amplification of roundoff errors when the solution is
reconstructed. The complexity of the algorithm is also analyzed, showing that
the cost of computing the spectral density function at a point grows less
rapidly than any fractional inverse power of the desired accuracy. A WKB
analysis is used to prove that the spectral density function is real analytic.
Using this new algorithm, we highlight key properties of the partial
differential equation and its solution that have strong implications on the
optimal choice of discretization method in large-scale plasma physics
computations.
