We present formulas for accurate numerical conversion between functions represented by multiwavelets and their multipole/local expansions with respect to the kernel of the form, \(e^{\lambda r}/r\). The conversion is essential for the application of fast multipole methods for functions represented by multiwavelets. The corresponding separated kernels exhibit near-singular behaviors at large \(\lambda\). Moreover, a multiwavelet basis function oscillates more wildly as its degree increases. These characteristics in combination render any brute-force approach based on numerical quadratures impractical. Our approach utilizes the series expansions of the modified spherical Bessel functions and the Cartesian expansions of solid harmonics so that the multipole-multiwavelet conversion matrix can be evaluated like a special function. The result is a quadrature-free, fast, reliable, and machine precision accurate scheme to compute the conversion matrix with predictable sparsity patterns.