For a specified hyperexponential function $F\left(x\,y\right)$ of x and y, the Zeilberger(F, x, y, Dx) calling sequence constructs for $F\left(x\,y\right)$ a Z-pair $L,G$ that consists of a linear differential operator with coefficients that are polynomials of x over the complex number field

and a hyperexponential function $G\left(x\,y\right)$ of x and y such that

Dx and Dy are the differential operators with respect to x, and y, respectively, defined by $\mathrm{Dx}\left(F\left(x\,y\right)\right)=\frac{\partial }{\partial x}F\left(x\,y\right)$, and $\mathrm{Dy}\left(F\left(x\,y\right)\right)=\frac{\partial }{\partial y}F\left(x\,y\right)$.

By assigning values to the global variables _MINORDER and _MAXORDER, the algorithm is restricted to finding a Z-pair $L,G$ for $F\left(x\,y\right)$ such that the order of L is between _MINORDER and _MAXORDER.

The algorithm has two implementations. The default implementation uses a variant of Gosper's algorithm, and another one is based on the universal denominators. With the 'gosper_free' option, Gosper-free implementation is used.

The output from the Zeilberger command is a list of two elements $\left[L\,G\right]$ representing the computed Z-pair $L,G$.

$\mathrm{with}\left(\mathrm{DEtools}\right)\:$

$F≔\exp \left(-\frac{x^2}{y^2}-y^2\right)$

$F≔{\ⅇ}^{-\frac{x^2}{y^2}-y^2}$

$\mathrm{Zpair}≔\mathrm{Zeilberger}\left(F\,x\,y\,\mathrm{Dx}\right)\:$

$L≔\mathrm{Zpair}\left[1\right]$

$L≔{\mathrm{Dx}}^2-4$

$G≔\mathrm{Zpair}\left[2\right]$

$G≔\frac{2{\ⅇ}^{-\frac{y^4+x^2}{y^2}}}{y}$

Almkvist, G, and Zeilberger, D. "The method of differentiating under the integral sign." Journal of Symbolic Computation. Vol. 10. (1990): 571-591.