I am skeptical about a possible closed form but, looking here, you will notice that there are quite good approximations
$$\left(\text{erf}(x)\right)^2\approx 1-e^{-a x^2}\qquad \text{with}\qquad a=(1+\pi )^{2/3} \log ^2(2)$$
So, you can approximate
$$I_n=\int_{-1}^1 \left(\text{erf}(x)\right)^{2n}\,dx$$by
$$J_n=\int_{-1}^1 \left(1-e^{-a x^2}\right)^{n}\,dx$$ for which the binomial expansion would be required (easy).
This would give you things like
$$J_1=2-\frac{\sqrt{\pi } \text{erf}\left(\sqrt{a}\right)}{\sqrt{a}}$$
$$J_2=2-\frac{2 \sqrt{\pi }
\text{erf}\left(\sqrt{a}\right)}{\sqrt{a}}+\frac{\sqrt{\frac{\pi }{2}}
\text{erf}\left(\sqrt{2a}\right)}{\sqrt{a}}$$
$$J_3=2-\frac{3 \sqrt{\pi } \text{erf}\left(\sqrt{a}\right)}{\sqrt{a}}+\frac{3
\sqrt{\frac{\pi }{2}} \text{erf}\left(\sqrt{2a}
\right)}{\sqrt{a}}-\frac{\sqrt{\frac{\pi }{3}} \text{erf}\left(\sqrt{3a}
\right)}{\sqrt{a}}$$
$$\color{blue}{J_n=2+\sqrt{\frac \pi a}\,\sum_{k=1}^n (-1)^k \frac{\binom{n}{k}}{\sqrt{k}}\,\text{erf}\left(\sqrt{ak} \right) }\tag 1$$ I produce below a short table for comparison
$$\left(
\begin{array}{ccc}
n & \text{approximation} & \text{exact} \\
1 & 0.591506 & 0.596751 \\
2 & 0.279674 & 0.283168 \\
3 & 0.151067 & 0.153256 \\
4 & 0.0870954 & 0.0884650 \\
5 & 0.0522216 & 0.0530855 \\
6 & 0.0321485 & 0.0326982 \\
7 & 0.0201718 & 0.0205243 \\
8 & 0.0128409 & 0.0130686 \\
9 & 0.00826756 & 0.00841548 \\
10 & 0.00537202 & 0.00546863
\end{array}
\right)$$
Update
After my answer to your next question, I reused for this problem my approach with the same Padé approximants and obtained as approximations
$$I_n=\frac 2{2n+1} \left(\frac{4}{\pi }\right)^n\,\, _2F_1\left(2 n,\frac{2n+1}{2};\frac{2n+3}{2};-\frac{1}{3}\right)\tag 2$$
$$I_n=\frac 2{2n+1} \left(\frac{4}{\pi }\right)^n\,F_1\left(\frac{2n+1}{2};-2 n,2 n;\frac{2n+3}{2};\frac{1}{30},-\frac{3}{10}\right)\tag 3$$
