La antiderivada se puede encontrar con la mano. Integrar por partes
$$\newcommand{\fer}{\operatorname{fer}}
\mathrm{d}u=\pecado 2\theta\,\mathrm{d}\theta,\quad v=\erf(\cos \theta)\erf(\sin \theta)
$$
da
$$
-\frac12\cos 2\theta\erf(\sin \theta)\erf(\cos \theta)+
\frac1{\sqrt\pi}\int\cos(2\theta)\left[e^{-\sin^2 \theta}\cos \theta\erf(\cos \theta)-e^{-\cos^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta
$$
Por lo tanto estamos reducidos a la computación
$$
\int\cos(2\theta)\left[e^{-\sin^2 \theta}\cos \theta\erf(\cos \theta)-e^{-\cos^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta\\
=e^{-1}\int\cos(2\theta)\left[e^{\cos^2 \theta}\cos \theta\erf(\cos \theta)-e^{\sin^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta
$$
Pero tenga en cuenta que
\begin{align*}
&\frac{\mathrm{d}}{\mathrm{d}\theta}\left[e^{\cos^2 \theta}\sin \theta\erf(\cos \theta)\right]\\
&=(-2\sin \theta\cos \theta)e^{\cos^2 \theta}\sin \theta\erf(\cos \theta)
+e^{\cos^2 \theta}\cos \theta\erf(\cos \theta)
+e^{\cos^2 \theta}\sin \theta\frac2{\sqrt\pi}e^{-\cos^2 \theta}(-\sin \theta)\\
&=(-2\sin^2 \theta+1)e^{\cos^2 \theta}\cos \theta\erf(\cos \theta)
-\frac{2\sin^2 \theta}{\sqrt\pi}\\
&=\cos (2 \theta)e^{\cos^2 \theta}\cos \theta\erf(\cos \theta)
-\frac{2\sin^2 \theta}{\sqrt\pi}
\end{align*}
y la ecuación similar con $\cos \theta$ e $\sin \theta$intercambiados.
Por lo que poner a estos dos juntos,
\begin{align*}
&\int\cos(2\theta)\left[e^{\cos^2 \theta}\cos \theta\erf(\cos \theta)-e^{\sin^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta\\
&=\left( e^{\cos^2 \theta} \erf(\cos \theta) \sin \theta - e^{\sin^2 \theta} \erf(\sin \theta) \cos \theta\right)+\int\frac2{\sqrt\pi}(\sin^2 \theta+\cos^2 \theta)\,\mathrm{d}\theta
\end{align*}
y, por supuesto, usted sabe cómo hacer el final integral.
El resultado final es
\begin{align*}
&\int\sin(2\theta)\erf(\sin\theta)\erf(\cos\theta)\,\mathrm{d}\theta\\
&=\frac2{e\pi}\theta-\frac12\cos 2\theta\erf(\sin \theta)\erf(\cos \theta)\\
&\quad+\frac1{e\sqrt\pi}
\left( e^{\cos^2 \theta} \erf(\cos \theta) \sin \theta - e^{\sin^2 \theta} \erf(\sin \theta) \cos \theta\right)
\end{align*}
Para $I$, este enfoque no funciona, y creo que no hay un closd forma de solución a la integral definida. Sin embargo, podemos expresar la serie de MacLaurin
$$
\tanh z=\sum_{k\geq 0} 2\frac{(-1)^k}{\pi^{2k+2}}\left(1-\frac1{4^{k+1}}\right)\zeta(2k+2)z^{2k+1}
$$
así
$$
I=4\sum_{k,k'\geq 0}(1-4^{-k-1})(1-4^{-k'-1})\frac{(-1)^{k+k'}\zeta(2k+2)\zeta(2k'+2)\Gamma(k+\frac12)\Gamma(k'+\frac12)}{\Gamma(k+k'+1)\pi^{2(k+k')+4}}
$$
que realmente no simplificar la tarea.