► The paper. The polynomials Ba_{n}^{(p)}(x) are studied in the paper "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/diffusionwave problems" by S. Bravo Yuste and E. Abad J. Phys. A: Math. Theor. 44 (2011) 075203. Ba_{n}^{(p)}(x) is a polynomial of degree 2n with an independent coefficient equal to 1.
► The name. To call Badajoz polynomials to these functions makes sense because Ba is the abbreviation of the name of the city and province of Badajoz ˗˗ the city where E. Abad and I live. To call them Bessel approximation polynomials or even BravoAbad polynomials, could be suitable, but I think that "Badajoz polynomials" is a nicer name.
► Two key properties of the Badajoz polynomials
►The normalized Bessel function. The relation between Bessel functions J_{p}(x) and Badajoz polynomials Ba_{n}^{(p)}(x) can be better appreciated by means of the normalized Bessel function :
↔
z_{p} being the first zero of J_{p}(x).


Bessel functions of the first kind for p=0,1,2 
Normalized Bessel functions for p=0,1,2 
►Ba_{n}^{(p)}(x) is an approximation to the normalized Bessel function. The polynomials Ba_{n}^{(p)}(x) are polynomial approximations to the normalized Bessel function , and consequently , also polynomial approximations to the Bessel function J_{p}(x):
↔
This approximation improves when n increases:
↔
● The Badajoz function J_{p,n}(x) defined above is then a polynomial approximation to the Bessel function of the first kind J_{p}(x). This Badajoz function can be recursively calculated by means of this iterative formula:
with
,
,
,
►The generating integral operators and . Families of increasingly better approximation functions to can be obtained by means of the repeated application of the integral operator over a seed function f_{0}(x). This operator is
where
When the seed function is equal to one, f_{0}(x)=1, we get the family of Badajoz polynomials:
For other seed functions f_{0}(x), we get other families of approximations. For example, for f_{0}(x)=1x we get the polynomials Be_{n}^{(p)}(x) (see S. Bravo Yuste and E. Abad J. Phys. A: Math. Theor. 44 (2011) 075203).
► Integral operators and with Mathematica:
f0[x_]:=1 
Λ[0,p_]:=f0[x] 
Λ[n_,p_]:=zp^2*Integrate[(1/u^(2p+1))*Integrate[x^(2p+1)*Λ[n1,p],{x,0,u},Assumptions>Re[p]>1&&Re[u]>0&&Re[zp]>0],{u, xx,1},Assumptions>Element[xx,Reals]&&Re[xx]>=0&&Re[zp]>0]/.xx>x 
Λhat[n_,p_]:=(Lan=Λ[n,p];Lad=Lan/.x>0;Collect[Lan/Lad//Simplify,x]) 
These instructions generate Ba_{n}^{(p)}(x). For example, Λhat[2,p] generates Ba_{2}^{(p)}(x). One can get other families of approximations using other definitions for f0[x_]. For example, f0[x_]:=1x leads to the Be_{n}^{(p)}(x) family.
► Polynomials Ba_{n}^{(p)}(x) with Mathematica. An easy way to generate Ba_{n}^{(p)}(x) is by means of the following Mathematica instructions:
bag[x_,p_,0]:=1 ; 
bag[x_,p_,m_]:=bag[x,p,m]=(1)^m*p!*x^(2m)/(2^(2m) m! (m+p)!)Sum[(1)^k*p!/(2^(2k) k! (k+p)!)*bag[x,p,mk],{k,1,m}] 
Ba[x_,p_,n_]:=bag[x,p,n]/bag[0,p,n] 
The Mathematica function Ba[x_,p_,n_] provides the polynomial Ba_{n}^{(p)}(x).
►The first 20 polynomials Ba_{n}^{(0)}(x):
► Plot of Ba_{n}^{(0)}(x). Plot of the normalized Bessel function (dashed line) and the first 21 polynomials Ba_{n}^{(0)}(x) with n=0,1,2,...20, (solid lines):
► Plot of J_{0,n}(x). Plot of the Bessel function J_{0}(x) (dashed line) and the first 21 Badajoz functions J_{0,n}(x) with n=0,1,2,...20, (solid lines):
►The first 11 polynomials Ba_{n}^{(3/2)}(x):
► Plot of Ba_{n}^{(3/2)}(x). Plot of the normalized Bessel function (dashed line) and the first 11 polynomials Ba_{n}^{(3/2)}(x) with n=0,1,2,...10, (solid lines):
► Plot of J_{3/2,n}(x). Plot of the Bessel function J_{3/2}(x) (dashed line) and the first 21 Badajoz functions J_{3/2,n}(x) with n=0,1,2,...20, (solid lines):
Normal diffusion 
Subdiffusion 

ddimensional equation  
Boundary condition 
c(R,t)=0 

Initial condition  c(r,0)=c_{0}  
Solution as superposition of modes  
Temporal evolution of the modes 


Spatial form of the modes  
Plot of the first mode (j=1) for the 3D problem  
Plot of the second mode (j=2) for the 3D problem 