Solution (source code)

= Solution

Fix a compact $K\subset\Omega$ and choose $\delta>0$ such that a closed $3\delta$-neighbourhood of $K$ is compactly contained in $\Omega$. The hypothesis gives a common bound $M$ on that neighbourhood. Choose a smooth radial <mollifier> $\rho_\delta$ supported in $B_\delta(0)$ with integral one. The <mean value property for harmonic functions>, integrated against its radial weights, gives
$$
u_j(x)=\int\rho_\delta(x-y)u_j(y)\,dy
$$
in a neighbourhood of $K$. Differentiating this fixed-radius convolution gives the uniform interior estimate
$$
|\nabla u_j(x)|\leq M\|\nabla\rho_\delta\|_1\leq C_\rho M/\delta.
$$
Thus the sequence is equicontinuous there; taking its pointwise limit gives the same local modulus of continuity for $u$. A finite small-ball cover of $K$ now proves <uniform convergence> directly: choose the cover so both $u_j$ and $u$ vary by less than $\varepsilon$ inside each ball, then use convergence at its finitely many centres to make $|u_j-u|<\varepsilon$ there. The <triangle inequality> gives $\sup_K|u_j-u|<3\varepsilon$ for large $j$.

Finally, pass to the limit in the reproducing convolution identity on smaller compact sets. It yields $u=u*\rho_\delta$, so $u$ is smooth locally. Also, for each compactly supported smooth <test function> $\phi$, uniform convergence on its support gives
$$
\int_\Omega u\Delta\phi=\lim_j\int_\Omega u_j\Delta\phi=0.
$$
Hence its smooth Laplacian is zero. We have proved \b[$u_j\to u$ uniformly on every compact subset and $u$ is harmonic]. Using the continuous representatives in the pointwise hypothesis makes this conclusion pointwise, rather than only almost everywhere. This is <compact convergence of locally bounded harmonic functions>.