I find a possible way to solve it, assuming some properties of polynomials which I will state in the end: Write $p_j(z)=a_j+\sqrt{-1}b_j$ for $1\leq j\leq n$. Equating the real and imaginary parts, we see that solving $(1)$ is equivalent to solving $$\tag{2} \frac{\partial a_j}{\partial x_i}+\frac{\partial b_j}{\partial y_i} +\frac{\partial a_i}{\partial x_j}+\frac{\partial b_i}{\partial y_j}=2h_{i\overline{j}},\\ \frac{\partial b_j}{\partial x_i}-\frac{\partial a_j}{\partial y_i} -\frac{\partial b_i}{\partial x_j}+\frac{\partial a_i}{\partial y_j}=0,\\ \frac{\partial a_j}{\partial x_i}+\frac{\partial b_j}{\partial y_i} -\frac{\partial a_i}{\partial x_j}-\frac{\partial b_i}{\partial y_j}=0 $$ where $1\leq i, j\leq n$. Note that $(2)$ is a system of linear PDEs with constant coefficients. Let \begin{equation} P_{ij}^{k}(X_1,X_2,...,X_{2n-1},X_{2n})=\left\{ \begin{array}{ll} X_{2j-1}, & \hbox{ if $k=2i-1$;} \\ X_{2j}, & \hbox{ if $k=2i$;}\\ X_{2i-1}, & \hbox{ if $k=2j-1$;} \\ X_{2i}, & \hbox{ if $k=2j$;}\\ 0, & \hbox{ otherwise.} \end{array} \right. \end{equation} and $$(u_1,u_2,...,u_{2n-1},u_{2n})=(a_1,b_1,...,a_n,b_n).$$ Then the system of equations in the first line of $(2)$ can be written as $$\sum_{1\leq k\leq l\leq n}P_{ij}^{k}\Big(\frac{\partial}{\partial x_1},\frac{\partial}{\partial y_1},...,\frac{\partial}{\partial x_n},\frac{\partial}{\partial y_n}\Big)u_k=0$$ where $1\leq i\leq j\leq n$.

I guess that if $\{Q_{i,j}\}_{1\leq i\leq j\leq n}$ are polynomials in $X_1$, $X_2$,..., $X_{2n-1}$, $X_{2n}$ such that $$ \sum_{1\leq i\leq j\leq n}Q_{i,j}P_{ij}^k=0\mbox{ for all }1\leq k\leq 2n, $$ then we must $Q_{i,j}=0$ for all $1\leq i\leq j\leq n$, which I have posted it as a question here. If this is true, we can then apply Theorem 7.6.13 in the book of Hormander [ An introduction to complex analysis in several variables, Third edition, North-Holland, Amsterdam, 1990] to conclude that the system $(2)$ is solvable.

I find a possible way to solve it, assuming some properties of polynomials which I will state in the end: Write $p_j(z)=a_j+\sqrt{-1}b_j$ for $1\leq j\leq n$. Equating the real and imaginary parts, we see that solving $(1)$ is equivalent to solving $$\tag{2} \frac{\partial a_j}{\partial x_i}+\frac{\partial b_j}{\partial y_i} +\frac{\partial a_i}{\partial x_j}+\frac{\partial b_i}{\partial y_j}=2h_{i\overline{j}},\\ \frac{\partial b_j}{\partial x_i}-\frac{\partial a_j}{\partial y_i} -\frac{\partial b_i}{\partial x_j}+\frac{\partial a_i}{\partial y_j}=0,\\ \frac{\partial a_j}{\partial x_i}+\frac{\partial b_j}{\partial y_i} -\frac{\partial a_i}{\partial x_j}-\frac{\partial b_i}{\partial y_j}=0 $$ where $1\leq i, j\leq n$. Note that $(2)$ is a system of linear PDEs with constant coefficients. Let \begin{equation} P_{ij}^{k}(X_1,X_2,...,X_{2n-1},X_{2n})=\left\{ \begin{array}{ll} X_{2j-1}, & \hbox{ if $k=2i-1$;} \\ X_{2j}, & \hbox{ if $k=2i$;}\\ X_{2i-1}, & \hbox{ if $k=2j-1$;} \\ X_{2i}, & \hbox{ if $k=2j$;}\\ 0, & \hbox{ otherwise.} \end{array} \right. \end{equation} and $$(u_1,u_2,...,u_{2n-1},u_{2n})=(a_1,b_1,...,a_n,b_n).$$ Then the system of equations in the first line of $(2)$ can be written as $$\sum_{1\leq k\leq l\leq n}P_{ij}^{k}\Big(\frac{\partial}{\partial x_1},\frac{\partial}{\partial y_1},...,\frac{\partial}{\partial x_n},\frac{\partial}{\partial y_n}\Big)u_k=0$$ where $1\leq i\leq j\leq n$.

I guess that if $\{Q_{i,j}\}_{1\leq i\leq j\leq n}$ are polynomials in $X_1$, $X_2$,..., $X_{2n-1}$, $X_{2n}$ such that $$ \sum_{1\leq i\leq j\leq n}Q_{i,j}P_{ij}^k=0\mbox{ for all }1\leq k\leq 2n, $$ then we must $Q_{i,j}=0$ for all $1\leq i\leq j\leq n$, which I have posted it as a question here. If this is true, we can then apply Theorem 7.6.13 in the book of Hormander [ An introduction to complex analysis in several variables, Third edition, North-Holland, Amsterdam, 1990] to conclude that the system $(2)$ is solvable.