This question is closely related to the following MO question Characterizing the real analytic Eisenstein series

Let $\mathfrak{h}=\{z=x+iy\in\mathbf{C}\}$ be the Poincare upper half plane endowed with its Poincare metric. Let $w\in\mathbf{Z}$ be a weight and define the "weight $w$ Laplacian" on $\mathfrak{h}$ by $$ \Delta[w]=-y^2(\partial_x^2+\partial_y^2)+i\cdot w\cdot y\cdot\partial_x. $$ Here $i=\sqrt{-1}$. In general, $\Delta[w]$ will commute with the $|_{w}$ right action of $GL_2^+(\mathbf{R})$ on Maass forms as is explained on page 130 of Bump's book on automorphic forms and representations. Consider the differential equation $$ (1)\;\;\;\;\;\;\;\; \Delta[w] f(y,s)=(s-w/2)(1-s+w/2)f(y), $$ where $f(y,s)$ is a function depending only on $y=Im(z)$ and $s\in\mathbf{C}$ and is such that $s\mapsto f(y,s)$ is meromorphic. Solving the corresponding linear ODE of order $2$ in $y$ we see that the solution space of (1) (for $s\neq \frac{1}{2}$) has the form $$ A(s)y^{s-w/2}+B(s)y^{1-s+w/2}. $$ where $A(s)$ and $B(s)$ are meromorphic functions in $s$.

Let $g\in\mathbf{Z}_{\geq 1}$ be a fixed integer and consider the symmetric space $\mathfrak{h}^g$. Let $\underline{z}=(z_j)_{j=1}^g\in\mathfrak{h}^g$ where $z_j=x_j+iy_j$. Let $\underline{w}=(w_1,w_2,\ldots,w_g)\in\mathbf{Z}^g$ be an integral weight vector. We let $D_j[w_j]$ be the weight $w_j$ Laplacian with respect to the variable $z_j$.

Consider now the following differential system of linear PDEs: $$ (2)\;\;\;\;\;\;\;\;\; D_j[w_j]f(\underline{y},s)=(s-w_j/2)(1-s+w_j/2)f(\underline{y},s)\;\;\;\;\;\;\;\;\; (1\leq j\leq g) $$ where $f(\underline{y},s)$ is a function depending only on $Im(\underline{z})=\underline{y}=(y_1,\ldots,y_g)$ and $s$ and which is such that $s\mapsto f(\underline{y},s)$ is meromorphic. Let $J\cup J'=\{1,2,\ldots,g\}$ with $J\cap J'=\emptyset$ be a partition of $\{1,2,\ldots,g\}$ in two blocks. A solution for $(2)$ is given by $$ (3)\;\;\;\;\;\;\;\;A_{I,I'}(s)\prod_{j\in J } y_j^{s-w_j/2}\prod_{j\in J'} y_j^{1-s+w_j/2}. $$ Here $A_{J,J'}(s)$ are meromorphic functions in $2$. Let us take the $\mathbf{C}$-vector space generated by all such solutions and call it $V$. Note that $V$ is a vector space of dimension $2^g$ over the ring of meromorphic functions in one variable.

Q: Are all the solutions of the differential system $(2)$ elements of $V$?

Remarks:

(a) Since $f(\underline{y},s)$ depends only on $\underline{y}$ and not on the variables $\underline{x}$, the action of $D_{j}[w_j]$ on $f(\underline{y})$ simplifies. Therefore the system (2) becomes a system of linear ODES of order $2$ in the variables $y_i$'s. But I preferred to keep $D_{j}[w_j]$ as I thought of it orginally, since I find it better motivated.

(b) So my heuristic reasoning which makes me think that the question has an affirmative answer is the following: Let us fix the index $j=1$ and let $f(\underline{y},s)$ be a solution to
$$ D_1[w_1]f(\underline{y},s)=(s-w_j/2)(1-s+w_j/2)f(\underline{y},s). $$ Then we know from ODE theory that $$ f(\underline{y},s)=A_1(s,y_2,\ldots,y_g)y_1^{s-w_j/2}+B_1(s,y_2,\ldots,y_g)y_1^{1-s+w_j/2} $$ For some functions $A_1(s,y_2,\ldots,y_g)$ and $B_1(s,y_2,\ldots,y_g)$. And then we can keep repeating this procedure until we sieve out all the variables $y_i$'s. But I feel a bit uneasy with this kind of reasoning because of a lack of personal experience with PDE systems...