Computing $_2F_2(a,a,a+1,a+1,z)$ (hypergeometric function) Trying to implement the derivative of the gamma incomplete function, I encountered the hypergeometric function $_2F_2(a,a,a+1,a+1; z=-x)$, where $x$ would always be a positive real (and thus $z$ a negative real), and $a$ could theoretically be any complex, but I only need positives reals too. The serie collapses to the simple :
$$a^2 \sum_{n = 0}^{\infty} \left(\frac{1}{a+n}\right)^2 \frac{z^n}{n!}.$$
Is there something I could do analytically to transform it (we might assume $a>0, z <0$ both reals) into something easier to compute ? Is there Litterature about this serie ?
 A: This is a partial answer for integer $a\in\mathbb N$, where the given hypergeometric function becomes a polynomial:
Defining the exponential generating function
$$
G_m(x)=\sum_{n=0}^m H_n \frac{x^n}{n!}\tag{1}
$$
of the harmonic numbers
$$
H_n = \sum_{k=1}^n \frac 1 k,\tag{2}
$$
which fulfills
$$
G_\infty(x)=e^x\big(\log(x) - \mathrm{Ei}(-x) + \gamma\big)\tag{3}
$$
the polynomial is given by
\begin{align}
\frac{{}_2F_2(\cdots;-x)}{a^2}&=
(1-\delta_{a,1})\frac{e^{-x}}{x^2}\\
&+\frac{(a-1)!}{x^a}
\left[
\left(\frac{\Gamma(a-2,x)}{(a-3)!}-1\right)H_{a-1} 
+e^{-x}\big(G_\infty(x)-G_{a-3}(x)\big)
\right]. \tag{4}
\end{align}
Here, $\mathrm{Ei}(-x)=-\Gamma(0,x)$ is the exponential integral function and $\Gamma(a,x)$ is the incomplete gamma function, which also reduces to a polynomial for integer $a$.
For $a=\{1,2\}$ this expression gives the terms in the comments of @Robert.
A: If numerical integration is acceptable for your application then you can use DLMF 16.5.2 to recast your hypergeometric function into a double integral on the unit square as
$$
{_2F_2}\left({a,a\atop a+1,a+1};z\right)=a^2\int_0^1\int_0^1t^{a-1}u^{a-1}e^{ztu}\,\mathrm dt\mathrm du.
$$
This form has the advantage of requiring no special functions in the integrand and, according to the DLMF, only requires $\Re a>0$.
If you do not have access to numerical integration but can use pseudo-random numbers, we can also obtain a Monte Carlo estimator for this hypergeometric function via
$$
{_2F_2}\left({a,a\atop a+1,a+1};z\right)=\mathsf Ee^{zTU},\quad T,U\sim\operatorname{Beta}(a,1).
$$
However, this will require $a$ to also be real.
A: To numerically evaluate this function, the quarter plane of possible $a$ and $x$ values will have to split into different regions. Near the origin, a truncation of the Maclaurin series can be effective. Asymptotic expansion can be effective for large $x$ and Mathematica gives me the following as an example:
$$e^{-x+O\left(\left(\frac{1}{x}\right)^2\right)}
   \left(\frac{a^2}{x^2}+O\left(\left(\frac{1}{x}\right)^3\right)\right)+x^{-a} \left(a
   \Gamma (a+1) (\log (x)-\psi ^{(0)}(a))+O\left(\left(\frac{1}{x}\right)^2\right)\right)$$
It's possible to generate more terms. Careful analysis must be done to verify in what regions these approximations can be useful.
(12) from this paper https://www.researchgate.net/publication/348049405_The_Computation_of_a_2F2_Hypergeometric_Function gives us an integral representation to play with too:
$$\, _2F_2(a,a;a+1,a+1;-x)=\frac{a^2}{x^a} \int_{0}^{1} \frac{\gamma(a,xt)}{t} dt $$
Some clever quadrature technique applied to this integral could be useful for an intermediate region.
Another step missing is to investigate what happens when $a$ is large. Starting with the integral representation is often fruitful for dealing with large parameters (although I have not yet found anything interesting with some quick Mathematica experiments...)
