On a novel iterative method to compute polynomial approximations to Bessel functions of the first kind and its connection to the solution of fractional diffusion/diffusion-wave problems

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 useful not only 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.