Solution (source code)

= Solution

For every compactly supported smooth <test function> $\phi$, integration by parts in the definition of the <weak derivative> gives
$$
\int f_j\partial_\ell\phi=-\int(\partial_\ell f_j)\phi.
$$
Both $\phi$ and $\partial_\ell\phi$ belong to $L^2(\mathbb R^n)$, so the assumed weak limits imply
$$
\int g\partial_\ell\phi=-\int h_\ell\phi.
$$
This is exactly the distributional characterization of $\partial_\ell g=h_\ell$. Since $g,h_1,\ldots,h_n$ are all in $L^2$, the defining condition of the <H1 space> is satisfied. Thus
$$
\boxed{g\in H^1(\mathbb R^n),\qquad\partial_\ell g=h_\ell\quad(1\leq\ell\leq n).}
$$
For complex functions the same argument uses conjugated test functions in the $L^2$ pairing. No prior identification of an $H^1$ weak limit is needed.