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.
