  We prove that the homogeneously polyanalytic functions of total order $m$,
defined by the system of equations $\overline{D}^{(k_1,\ldots,k_n)} f=0$ with
$k_1+\cdots+k_n=m$, can be written as polynomials of total degree $<m$ in
variables $\overline{z_1},\ldots,\overline{z_n}$, with some analytic
coefficients. We establish a weighted mean value property for such functions,
using a reproducing property of Jacobi polynomials. After that, we give a
general recipe to transform a reproducing kernel by a weighted change of
variables. Applying these tools, we compute the reproducing kernel of the
Bergman space of homogeneously polyanalytic functions on the unit ball in
$\mathbb{C}^n$ and on the Siegel domain. For the one-dimensional case,
analogous results were obtained by Koshelev (1977), Pessoa (2014), Hachadi and
Youssfi (2019).
