Where's my mistake in my attempt at showing that the squared sum of normally distributed variables is a $\chi^2$ distribution?

In your second equality you don't need to do a change of variables. You end up doing another change of variables when switching to $u$, which is the only substitution you actually need to do. \begin{align} P(X_1^2 \in [a^2, b^2]) &= P(X_1 \in [a,b] \cup [-b, -a]) \\ &= \int_{[a,b] \cup [-b, -a]} f_{X_1}(x) \, dx \\ &= 2 \int_a^b f_{X_1}(x) \, dx & \text{symmetry} \\ &= 2 \int_{a^2}^{b^2} \frac{1}{2\sqrt{u}} f_{X_1}(\sqrt{u}) \, du & u=x^2 \\ &= \int_{a^2}^{b^2} \frac{1}{\sqrt{2\pi}} u^{-1/2} e^{-u/2} \, du. \end{align}

Regarding upgrading "agreement on integrals" to "agreement on Borel sets," see Carathéodory's Extension Theorem.

$$\begin{aligned}P(X^2_1\in[a,b])&=P(X_1\in [-\sqrt{b},-\sqrt{a}]\cup [\sqrt{a},\sqrt{b}])=\\ &=P(X_1\in[-\sqrt{b},-\sqrt{a}])+P(X_1\in [\sqrt{a},\sqrt{b}])=\\ &=2P(X_1\in [\sqrt{a},\sqrt{b}])=\\ &=2(\Phi(\sqrt{b})-\Phi(\sqrt{a}))\end{aligned}$$ At this point, find the cdf to obtain the pdf by differentiation. Now let $a\to -\infty$ we get $$F_{X_1^2}(b)=2\Phi(\sqrt{b})-1\implies f_{X_1^2}(b)=\frac{1}{\sqrt{b}}\phi(\sqrt{b})=\frac{1}{\sqrt{2\pi b}}e^{-\frac{b}{2}}$$ This is the $\chi^2$ pdf with $n=1$. As suggested by @Bey, now use CFs. The CF of a $Y_1\sim \chi^2$ is $$\begin{aligned}E[e^{i\xi Y_1}]=(1-2i\xi)^{-1/2} \end{aligned}$$ So if we have the sum $Y=\sum_{k\leq n}Y_k$ where $Y_k$ are IID $\chi^2$ with $n=1$ we get $$E[e^{i\xi Y}]=E[e^{i\xi Y_1}]^n=(1-2i\xi)^{-n/2}$$ this means that $Y \sim \chi^2_n$.