If $\Delta u = f$ in $B_1$ where $f$ is a continuous function:
(a) then there exists a dimensional constant $C_n > 0$ such that $$\sup_{B_{\frac{1}{2}}} |\nabla u| \leq C_n \left( \sup_{B_1} |f| + \sup_{\partial B_1} |u| \right).$$ (b) If in addition, $\Delta u_n = f_n$ in $B_1$ and $|u_n| \leq C$ and $f_n \to 0$ uniformly then (up to a subsequence) $u_n$ converges to a harmonic function $h$.
For part (a) I was thinking of using the solution formula for the Dirichlet problem on the unit ball, that is:
\begin{align*} u(x) = \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \end{align*} Then proceeding by taking the gradient, which yields:
\begin{align*} \nabla u(x) &= \nabla \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \\ &= \left( \nabla \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \right) \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y + \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \nabla \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \\ &= -\frac{\nabla \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y + \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} \nabla \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \\ &= -\frac{2 \lvert x \rvert \nabla \lvert x \rvert}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y + \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} -n \lvert x - y \rvert^{-n - 1} \nabla \lvert x - y \rvert u(y) \mathrm{d}y \\ &= -\frac{2 x \nabla x}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y + \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} -n \lvert x - y \rvert^{-n - 2} (x - y) \nabla (x - y) u(y) \mathrm{d}y \\ &= -\frac{2 x}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y + \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} -n \lvert x - y \rvert^{-n - 2} (x - y) u(y) \mathrm{d}y \\ &= -\frac{2 x}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} u(y) \mathrm{d}y + \frac{1 - \lvert x \rvert^2}{\sigma_{n-1}} \int_{\partial B_1} -n \lvert x - y \rvert^{-n} \frac{x - y}{\lvert x - y \rvert^2} u(y) \mathrm{d}y \\ &= -\frac{1}{\sigma_{n-1}} \int_{\partial B_1} 2x \lvert x - y \rvert^{-n} u(y) \mathrm{d}y - \frac{1}{\sigma_{n-1}} \int_{\partial B_1} (1 - \lvert x \rvert^2) n \lvert x - y \rvert^{-n} \frac{x - y}{\lvert x - y \rvert^2} u(y) \mathrm{d}y \\ &= \frac{1}{\sigma_{n-1}} \int_{\partial B_1} \left( -2x - n \frac{x - y}{\lvert x - y \rvert^2} \right) \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \\ &= -\frac{1}{\sigma_{n-1}} \int_{\partial B_1} \left( 2x + n \frac{x - y}{\lvert x - y \rvert^2} \right) \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \\ \end{align*}
Hence we have
\begin{align*} \lvert \nabla u(x) \rvert &= \frac{1}{\sigma_{n-1}} \left\lvert \int_{\partial B_1} \left( 2x + n \frac{x - y}{\lvert x - y \rvert^2} \right) \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \right\rvert \\ &\leq \frac{1}{\sigma_{n-1}} \int_{\partial B_1} \left\lvert \left( 2x + n \frac{x - y}{\lvert x - y \rvert^2} \right) \lvert x - y \rvert^{-n} u(y) \mathrm{d}y \right\rvert \\ &\leq \frac{1}{\sigma_{n-1}} \int_{\partial B_1} \left\lvert 2x + n \frac{x - y}{\lvert x - y \rvert^2} \right\rvert \lvert x - y \rvert^{-n} \lvert u(y) \rvert \mathrm{d}y \\ &\leq \frac{1}{\sigma_{n-1}} \int_{\partial B_1} \left( \lvert 2x \rvert + \left\lvert n \frac{x - y}{\lvert x - y \rvert^2} \right\rvert \right) \lvert x - y \rvert^{-n} \lvert u(y) \rvert \mathrm{d}y \\ &\leq \frac{1}{\sigma_{n-1}} \int_{\partial B_1} \left( 2 \lvert x \rvert + \frac{n}{\lvert x - y \rvert} \right) \lvert x - y \rvert^{-n} \lvert u(y) \rvert \mathrm{d}y \\ &\leq \frac{2 \lvert x \rvert}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n} \lvert u(y) \rvert \mathrm{d}y + \frac{n}{\sigma_{n-1}} \int_{\partial B_1} \lvert x - y \rvert^{-n-1} \lvert u(y) \rvert \mathrm{d}y \end{align*}
But I am unsure how to proceed from here for part (a).
Now for part (b) I want to use the observation that there exists a constant $C<\infty$ such that for all $u$ harmonic and bounded by $M$ on any ball $B(x_0,2r)$ that $$|u(x) - u(x_0)| \leq \sup_{B_r(x_0)} |\nabla u| |x-x_0| \leq \dfrac{CM}{r}|x-x_0|$$ for all $x \in B(x_0,r)$. However, I am unsure how to tie that together with $\Delta u_n = f_n \to 0$ uniformly.