In https://en.wikipedia.org/wiki/Meijer_G-function it is stated:
With the modern definition, the majority of the established special functions can be represented in terms of the Meijer G-function. A notable property is the closure of the set of all G-functions not only under differentiation but also under indefinite integration. In combination with a functional equation that allows to liberate from a G-function $G(z)$ any factor $z^ρ$ that is a constant power of its argument $z$, the closure implies that whenever a function is expressible as a G-function of a constant multiple of some constant power of the function argument, $f(x) = G(cx^γ)$, the derivative and the antiderivative of this function are expressible so too.
So I wonder if there is a formula for $\int G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , cz^\gamma\right) \,dz$?
For $z>0$, $c\in\mathbb R \setminus\{0\}$, $\gamma>0$ we can substitute $u:=cz^\gamma$ to get $$\int G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , cz^\gamma\right) \,dz = \frac{1}{\gamma c^{1/\gamma}} \int G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , u\right) \, u^{\frac{1-\gamma}{\gamma}} du .$$
Then using $ z^\rho \, G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , z \right) = G_{p,q}^{m,n}\left(\begin{matrix}a_1+\rho , \ldots, a_n+\rho, a_{n+1}+\rho,\ldots,a_p+\rho\\ b_1+\rho,\ldots,b_m+\rho,b_{m+1}+\rho,\ldots,b_q+\rho\end{matrix} , z\right) $ with $\rho=\frac{1-\gamma}{\gamma}$ and $\int G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , z\right) \,dz = G_{p+1,q+1}^{m,n+1}\left(\begin{matrix}1, a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,0,b_{m+1},\ldots,b_q\end{matrix} , z\right) $ and substituting back yields
$$ \int G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , cz^\gamma\right) \,dz = \frac{1}{\gamma c^{1/\gamma}} G_{p+1,q+1}^{m,n+1}\left(\begin{matrix}1,a_1+\frac{1}{b}, \ldots, a_n+\frac{1}{b}, a_{n+1}+\frac{1}{b},\ldots,a_p+\frac{1}{b}\\ b_1+\frac{1}{b},\ldots,b_m+\frac{1}{b},0,b_{m+1}+\frac{1}{b},\ldots,b_q+\frac{1}{b}\end{matrix} , cz^\gamma\right). $$
Now, in the case of $\gamma$ positive and uneven and $G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , 0\right)=0$ the integral kernel $z\mapsto G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , cz^\gamma\right)$ viewed as function over the reals is symmetric with respect to the origin. Therefore, in this case the formula extends to negative $z\in\mathbb R$ as well. For $\gamma$ positive and with even parity and again $G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , 0\right)=0$ the integral kernel as function over the reals is symmetric with respect to the y-Axis. Therefore, in this case (interpreting the right hand side as $0$ for $z=0$)
$$ \int G_{p,q}^{m,n}\left(\begin{matrix}a_1, \ldots, a_n, a_{n+1},\ldots,a_p\\ b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix} , cz^\gamma\right) \,dz = \frac{z/\sqrt{z^2}}{ba^{1/b}} G_{p+1,q+1}^{m,n+1}\left(\begin{matrix}1,a_1+\frac{1}{b}, \ldots, a_n+\frac{1}{b}, a_{n+1}+\frac{1}{b},\ldots,a_p+\frac{1}{b}\\ b_1+\frac{1}{b},\ldots,b_m+\frac{1}{b},0,b_{m+1}+\frac{1}{b},\ldots,b_q+\frac{1}{b}\end{matrix} , cz^\gamma\right). $$
Is it true, that in those two cases the formula is even valid for general $z\in\mathbb C\setminus\{0\}$? What if $\gamma >0$ is neither an even nor uneven number but, let's say $3/2$?
Motivation: This question is related to https://github.com/sympy/sympy/issues/25786