  We present an iterative method to obtain approximations to Bessel functions
of the first kind $J_p(x)$ ($p>-1$) via the repeated application of an integral
operator to an initial seed function $f_0(x)$. The class of seed functions
$f_0(x)$ leading to sets of increasingly accurate approximations $f_n(x)$ is
considerably large and includes any polynomial. When the operator is applied
once to a polynomial of degree $s$, it yields a polynomial of degree $s+2$, and
so the iteration of this operator generates sets of increasingly better
polynomial approximations of increasing degree. We focus on the set of
polynomial approximations generated from the seed function $f_0(x)=1$. This set
of polynomials is not only useful for the computation of $J_p(x)$, but also
from a physical point of view, as it describes the long-time decay modes of
certain fractional diffusion and diffusion-wave problems.
