[proofplan]
We define the Newton potential $\tilde{u} := \Phi * f$ and show it is bounded. Any bounded solution $u$ of $-\Delta u = f$ then differs from $\tilde{u}$ by a bounded harmonic function $w = u - \tilde{u}$, which must be constant by [Liouville's Theorem](/theorems/38). The boundedness of $\tilde{u}$ uses two ingredients: for $|x|$ large, the compact support of $f$ and the decay $|x - y|^{-(n-2)} \le (|x|/2)^{-(n-2)}$ give $\tilde{u}(x) \to 0$; for $|x|$ bounded, $\Phi \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$ (valid since $n \ge 3$) and $f \in L^\infty$ ensure finiteness.
[/proofplan]
[step:Define the Newton potential and verify it is bounded on $\mathbb{R}^n$]
Define
\begin{align*}
\tilde{u}: \mathbb{R}^n &\to \mathbb{R} \\
x &\mapsto \int_{\mathbb{R}^n} \Phi(x - y)\, f(y)\, d\mathcal{L}^n(y).
\end{align*}
Since $f \in C_c^2(\mathbb{R}^n)$, there exists $R > 0$ with $\operatorname{supp} f \subset B(0, R)$.
**Boundedness for $|x| \ge 2R$:** For $y \in B(0, R)$ and $|x| \ge 2R$, the triangle inequality gives $|x - y| \ge |x| - |y| \ge |x|/2$, hence $|x - y|^{-(n-2)} \le (|x|/2)^{-(n-2)} = 2^{n-2} |x|^{-(n-2)}$. Therefore
\begin{align*}
|\tilde{u}(x)| &\le \int_{B(0,R)} \frac{\Gamma(1 + n/2)}{n(n-2)\pi^{n/2}}\, \frac{|f(y)|}{|x - y|^{n-2}}\, d\mathcal{L}^n(y) \\
&\le \frac{2^{n-2}\, \Gamma(1 + n/2)}{n(n-2)\pi^{n/2}}\, \frac{\|f\|_{L^1(\mathbb{R}^n)}}{|x|^{n-2}},
\end{align*}
which tends to $0$ as $|x| \to \infty$ (using $n \ge 3$, so $n - 2 \ge 1$).
**Boundedness for $|x| \le 2R$:** The fundamental solution $\Phi$ satisfies $\Phi \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$ for $n \ge 3$ (as verified in [Theorem 5](/theorems/5)), so $|\tilde{u}(x)| \le \|\Phi\|_{L^1(B(0, 3R))}\, \|f\|_{L^\infty(\mathbb{R}^n)} < \infty$ for all $|x| \le 2R$.
Combining the two regions, $\tilde{u}$ is bounded on $\mathbb{R}^n$.
[guided]
The goal is to show that the Newton potential $\tilde{u} = \Phi * f$ is bounded, so that the difference $u - \tilde{u}$ is a bounded harmonic function and Liouville's theorem applies.
Define
\begin{align*}
\tilde{u}: \mathbb{R}^n &\to \mathbb{R} \\
x &\mapsto \int_{\mathbb{R}^n} \Phi(x - y)\, f(y)\, d\mathcal{L}^n(y).
\end{align*}
Since $f \in C_c^2(\mathbb{R}^n)$, there exists $R > 0$ with $\operatorname{supp} f \subset B(0, R)$, so the integral reduces to $\int_{B(0,R)}$.
We split $\mathbb{R}^n$ into two regions. For $|x| \ge 2R$: every $y \in B(0, R)$ satisfies $|x - y| \ge |x| - |y| \ge |x| - R \ge |x|/2$, so $\Phi(x - y) \le \frac{\Gamma(1 + n/2)}{n(n-2)\pi^{n/2}} \cdot \frac{2^{n-2}}{|x|^{n-2}}$. This gives
\begin{align*}
|\tilde{u}(x)| \le \frac{2^{n-2}\, \Gamma(1 + n/2)}{n(n-2)\pi^{n/2}}\, \frac{\|f\|_{L^1(\mathbb{R}^n)}}{|x|^{n-2}} \to 0 \quad \text{as } |x| \to \infty.
\end{align*}
The decay exponent $n - 2 \ge 1$ is where the hypothesis $n \ge 3$ is used. In dimension $n = 2$, the fundamental solution $\Phi(x) = -\frac{1}{2\pi} \log|x|$ does not decay at infinity, and $\tilde{u}$ need not be bounded.
For $|x| \le 2R$: the integration takes place over $y \in B(0, R)$, so $|x - y| \le 3R$. The singularity $|x - y|^{-(n-2)}$ is locally integrable in $\mathbb{R}^n$ for $n \ge 3$ (as $n - 2 < n$), which is established in [Theorem 5](/theorems/5). Therefore $|\tilde{u}(x)| \le \|\Phi\|_{L^1(B(0, 3R))}\, \|f\|_{L^\infty} < \infty$.
Combining, $\tilde{u}$ is bounded on $\mathbb{R}^n$.
[/guided]
[/step]
[step:Show the difference $w = u - \tilde{u}$ is bounded and harmonic, then apply Liouville's theorem]
Let $u \in C^2(\mathbb{R}^n)$ be a bounded solution of $-\Delta u = f$. By [Theorem 5](/theorems/5), $\tilde{u} = \Phi * f$ satisfies $-\Delta \tilde{u} = f$ in $\mathbb{R}^n$. Define $w := u - \tilde{u}$. Then
\begin{align*}
-\Delta w = -\Delta u - (-\Delta \tilde{u}) = f - f = 0 \quad \text{in } \mathbb{R}^n,
\end{align*}
so $w$ is harmonic on $\mathbb{R}^n$. Moreover, $w$ is bounded since $u$ is bounded by hypothesis and $\tilde{u}$ is bounded by the previous step. By [Liouville's Theorem](/theorems/38), $w$ is constant: $w \equiv C$ for some $C \in \mathbb{R}$. Therefore
\begin{align*}
u(x) = \tilde{u}(x) + C = \int_{\mathbb{R}^n} \Phi(x - y)\, f(y)\, d\mathcal{L}^n(y) + C,
\end{align*}
which is the desired representation.
[/step]