Following the comment by whuber, for problems like this that involve convolutions of IID random variables, it is generally simpler to work with the characteristic function than with the density function. Using the law of the unconscious statistician we can obtain the characteristic function for $X^2$ as follows:
$$\begin{align}
\phi_{X^2}(t)
&\equiv \mathbb{E}(\exp(itX^2)) \\[16pt]
&= \int \limits_{-\infty}^\infty \exp(it x^2) \cdot \text{N}(x|0, \sigma^2) \ dx \\[6pt]
&= \int \limits_{-\infty}^\infty \exp(it x^2) \cdot \frac{1}{\sqrt{2 \pi \sigma^2}} \cdot \exp \bigg( -\frac{1}{2 \sigma^2} \cdot x^2 \bigg) \ dx \\[6pt]
&= \int \limits_{-\infty}^\infty \frac{1}{\sqrt{2 \pi \sigma^2}} \cdot \exp \bigg( -\frac{1-2it \sigma^2}{2 \sigma^2} \cdot x^2 \bigg) \ dx \\[6pt]
&= \frac{1}{\sqrt{1-2it \sigma^2}} \int \limits_{-\infty}^\infty \frac{1}{\sqrt{2 \pi}} \sqrt{\frac{1-2it \sigma^2}{\sigma^2}} \cdot \exp \bigg( -\frac{1-2it \sigma^2}{2 \sigma^2} \cdot x^2 \bigg) \ dx \\[6pt]
&= \frac{1}{\sqrt{1-2it \sigma^2}} \int \limits_{-\infty}^\infty \text{N}\bigg( x \bigg| 0, \frac{\sigma^2}{1-2it \sigma^2} \bigg) \ dx \\[6pt]
&= \frac{1}{\sqrt{1-2it \sigma^2}}. \\[6pt]
\end{align}$$
(And of course, we have $\phi_{X^2}(t) = \phi_{Y^2}(t)$ in this case so the latter characteristic function is the same.) We then have:
$$\begin{align}
\phi_Z(t)
&\equiv \mathbb{E}(\exp(itZ)) \\[16pt]
&= \mathbb{E}(\exp(itX^2 + itY^2)) \\[16pt]
&= \mathbb{E}(\exp(itX^2)) \cdot \mathbb{E}(\exp(itY^2)) \\[12pt]
&= \frac{1}{\sqrt{1-2it \sigma^2}} \cdot \frac{1}{\sqrt{1-2it \sigma^2}} \\[6pt]
&= \frac{1}{1-2it \sigma^2}. \\[6pt]
\end{align}$$
This is the characteristic function for the scaled chi-squared distribution with two degrees-of-freedom. Using the fact that the characteristic function is a unique representative of the distribution, you then have:
$$Z \sim \sigma^2 \cdot \text{ChiSq}(\text{df} = 2).$$