Plancherel's theorem

Step 1 (a special case). First check the identity for a Gaussian $f(x) = e^{-\pi x^2}$. A direct computation (completing the square in the exponent) gives $\hat f(\xi) = e^{-\pi \xi^2}$, so both sides of $\int_{-\infty}^{\infty} |f(x)|^2\,dx = \int_{-\infty}^{\infty} |\hat f(\xi)|^2\,d\xi$ equal the same Gaussian integral $\int e^{-2\pi x^2}\,dx$.

Step 2 (extend by linearity and the convolution theorem). For nice (Schwartz) functions $f$, write $g = f * \tilde f$ where $\tilde f(x) = \overline{f(-x)}$; then $g(0) = \int |f(x)|^2\,dx$, and by the convolution theorem $\widehat{f * g} = \hat f \cdot \hat g$ together with Fourier inversion, $g(0) = \int \hat g(\xi)\,d\xi = \int |\hat f(\xi)|^2\,d\xi$. This proves the identity for all Schwartz functions, which include Gaussians and all smooth rapidly-decaying functions.

Step 3 (density argument). Schwartz functions are dense in $L^2(\mathbb{R})$: every $f \in L^2$ is a limit $f_n \to f$ of Schwartz functions in the $L^2$ norm. Since the Fourier transform is an isometry on this dense subspace by Step 2, it extends uniquely to a bounded operator on all of $L^2(\mathbb{R})$ that still satisfies $\int_{-\infty}^{\infty} |f(x)|^2\,dx = \int_{-\infty}^{\infty} |\hat f(\xi)|^2\,d\xi$ — this extension is what 'the Fourier transform of an $L^2$ function' means when the defining integral $\hat f(\xi) = \int_{-\infty}^{\infty} f(x)\, e^{-2\pi i x \xi}\,dx$ may not converge absolutely.