The paper focuses on an approximation of the first passage time probability density function of a Feller stochastic process by using cumulants and a Laguerre-Gamma polynomial approximation. The feasibility of the method relies on closed form formulae for cumulants and moments recovered from the Laplace transform of the probability density function and using the algebra of formal power series. To improve the approximation, sufficient conditions on cumulants are stated. The resulting procedure is made easier to implement by the symbolic calculus and a rational choice of the polynomial degree depending on skewness, kurtosis and hyperskewness. Some case-studies coming from neuronal and financial fields show the goodness of the approximation even for a low number of terms. Open problems are addressed at the end of the paper.