We examine the regularity of weak solutions to the semi-linear bi-Laplacian equation \[ \Delta^2u=f(x,u,Du) \quad \text{ in } \Omega, \] where $\Omega \subset \mathbb{R}^d$ is a bounded smooth domain, and the nonlinearity $f: \Omega \times \mathbb{R}\times \mathbb{R}^d \to\mathbb{R}$ satisfies the growth condition depending on $u$ and $Du$. We obtain that $u\in C^{2,\sigma}_{\operatorname{loc}}(\Omega)$ for $\sigma\in (0,1)$. That means, \[
\sum_{|\gamma|\leq 2}\|D^{\gamma}u\|_{C(\Omega')} + \sum_{|\gamma|=2} \sup_{\substack{x, y \in \Omega' \\ x \neq y}} \left\{\frac{\left|D^\gamma u(x) - D^\gamma u(y)\right|}{|x - y|^\sigma}\right\}<\infty, \] where $\Omega'\Subset \Omega$. Our strategy is to render this fourth-order problem as a system of two Poisson equations and explore the interplay between the integrability and smoothness available for each equation taken isolated. This is joint work with Edgard A. Pimentel and José Miguel Urbano.