This is a special case of Integration of radial functions?.

I don't really understand the solution given there.

Also, I only need the special case of $n=3$, i.e. how to prove $$\int_{\mathbb{R}^3}f(|x|)\,dx=4\pi\int_0^\infty f(r)r^2\,dr$$ when $f$ is not necessarily continuous.

The context will be Lebesgue integration rather than Riemann integration. The question actually warned that " we need to justify usage of Polar Coordinates since $f$ may not be continuous".

Thanks for any help. A reference will be acceptable too; the Stromberg reference given in the previous answer was too general and difficult to understand for me.

What I know is this theorem which works in Lebesgue Integral: Let $y=Tx$ be a nonsingular transformation of $\mathbb{R}^n$. If $\int_Ef(y)\,dy$ exists, then

$\int_E f(y)\,dy=|\det T|\int_{T^{-1}E}f(Tx)\,dx$

One of the standard machinery in measure theory is the monotone class theorem and its variants. The essence of this line of reasoning is that you may first establish your claim for simple and nice cases and then extend the result by appealing to the standard monotone-class argument.

Let us apply this idea to our problem. Let $X = [0,\infty)$ be the closed half-interval and $\mathcal{B}(X)$ be the Borel $\sigma$-algebra on $X$. Then the following claim is true:

Claim. Suppose that $\mu$ and $\nu$ are measures on $\mathcal{B}(X)$ such that $$\mu([a,b]) = \nu([a,b]) < \infty \quad : \ \forall 0 \leq a \leq b < \infty.$$ Then $\mu = \nu$ on $\mathcal{B}(X)$.

The proof is an easy application of the Dynkin's $\pi$-$\lambda$ theorem. So let me explain how this answers to our problem.

Step 1. Define two measures $\mu$ and $\nu$ on $X = [0,\infty)$ by

$$\mu(E) = \int_{\Bbb{R}^3} \mathbf{1}_E(|x|) \, dx, \qquad \nu(E) = \int_{0}^{\infty} 4\pi r^2 \mathbf{1}_E (r) \, dr.$$

By direct computation, $\mu$ and $\nu$ agree when $E$ is any bounded closed interval:

$$\mu([a,b]) = \operatorname{Leb}(x \in \Bbb{R}^3 : a \leq |x| \leq b) = \frac{4\pi}{3}(b^3 - a^3) = \nu([a,b]).$$

By the claim, $\mu = \nu$ on $\mathcal{B}(X)$. Together with the fact that $\nu$ is absolutely continuous w.r.t. the Lebesgue measure on $X$, it easily follows that $\mu = \nu$ on the Lebesgue $\sigma$-algebra as well.

Step 2. The previous step shows that the desired change-of-variables formula holds for indicator functions. Then the general case follows by approximating $f$ by simple functions.

Alternatively, we may utilize the layer-cake representation to give another proof. Assume WLOG that $f \geq 0$. If we define $E(t) = \{ r \in X : t < f(r) \}$, then by the Tonelli's theorem,

\begin{align*} \int_{\Bbb{R}^3} f(|x|) \, dx &= \int_{\Bbb{R}^3} \left( \int_{0}^{\infty} \mathbf{1}_{\{t < f(|x|)\}} \, dt \right) \, dx \\ &= \int_{0}^{\infty} \int_{\Bbb{R}^3} \mathbf{1}_{\{t < f(|x|)\}} \, dxdt = \int_{0}^{\infty} \mu(E(t)) \, dt \end{align*}

and likewise

\begin{align*} \int_{0}^{\infty} 4\pi r^2 f(r) \, dr &= \int_{0}^{\infty} 4\pi r^2 \left( \int_{0}^{\infty} \mathbf{1}_{\{t < f(r)\}} \, dt \right) \, dr \\ &= \int_{0}^{\infty} \int_{0}^{\infty} 4\pi r^2 \mathbf{1}_{\{t < f(r)\}} \, drdt =\int_{0}^{\infty} \nu(E(t)) \, dt. \end{align*}

So they are equal.