This talk discusses the use of numerical methods to approximate the probability density function of random variables governed by compound distributions. These random variables are useful in actuarial science to model the risk of a portfolio of contracts. In ruin theory, the probability of ultimate ruin within the compound Poisson ruin model is the survival function of a geometric compound distribution. The proposed method consists in a projection of the probability density function onto an orthogonal polynomial system. These polynomials are orthogonal with respect to a probability measure that belongs to Natural Exponential Families with Quadratic Variance Function. The polynomial approximation is compared to other numerical methods that recover the probability density function from the knowledge of the moments or the Laplace transform of the distribution. The polynomial method is then extended in a multidimensional setting, along with the probability density estimator derived from the approximation formula.