Let
$$
S(t)=\int_{0}^{t}\sin(x^{2})\,dx\quad\text{and}\quad C(t)=\int_{0}^{t}\sin(x^{2})\,dx\,.
$$
Using the identities
$$
\sin(x^{2})=-\frac{1}{2x}\frac{d}{dx}\cos(x^{2})\quad\text{and}\quad\cos(x^{2})=\frac{1}{2x}\frac{d}{dx}\sin(x^{2})\,,
$$
integrating by parts over $[1,t]$, and passing to the limit as $t\to\infty$, one shows that there exist
$$
\lim_{t\to\infty}S(t)=S_{\infty}\in\mathbb{R}\quad\text{and}\quad\lim_{t\to\infty}C(t)=C_{\infty}\in\mathbb{R}\,.
$$
Moreover, since
$$
S_{\infty}=\int_{0}^{\infty}\frac{\sin y}{2\sqrt y}\,dy=\sum_{k=0}^{\infty}(-1)^{k}a_{k}\quad\text{where}\quad a_{k}=\int_{k\pi}^{(k+1)\pi}\frac{|\sin y|}{2\sqrt y}\,dy
$$
and $0<a_{k+1}<a_{k}$, we infer that $S_{\infty}>0$. To evaluate $S_{\infty}$ with the Feynman method, let us introduce the mappings
$$
f(t)=\left(\int_{0}^{t}\sin(x^{2})\,dx\right)^{2}+\left(\int_{0}^{t}\cos(x^{2})\,dx\right)^{2},\quad
g(t)=\int_{0}^{1}\frac{\sin(t^{2}(1-x^{2}))}{1-x^{2}}\,dx\,.
$$
We can check that $f'(t)=g'(t)$ and $f(0)=g(0)$, hence $f(t)=g(t)$ for every $t\ge 0$. Since
$$
g(t)=\int_{0}^{t^{2}}\frac{\sin(2x+x^{2}/t^{2})}{2x+x^{2}/t^{2}}\,dx\to\int_{0}^{\infty}\frac{\sin(2x)}{2x}\,dx=\frac{\pi}{4}\quad\text{as }t\to\infty\,,
$$
we obtain that
$$
S_{\infty}^{2}+C_{\infty}^{2}=\frac{\pi}{4}\,.
$$
Proving that
$$
(*)\qquad S_{\infty}^{2}=C_{\infty}^{2}\,,
$$
we conclude that $S_{\infty}=\sqrt{\pi/8}$. To show $(*)$ we introduce the mappings
$$
F(t)=\left(\int_{0}^{t}\cos(x^{2})\,dx\right)^{2}-\left(\int_{0}^{t}\sin(x^{2})\,dx\right)^{2},\quad
G(t)=\int_{0}^{1}\frac{\sin(t^{2}(1+x^{2}))}{1+x^{2}}\,dx\,.
$$
One can check that $F'(t)=G'(t)$ and $F(0)=G(0)$, hence $F(t)=G(t)$ for every $t\ge 0$ and in particular
$$
C_{\infty}^{2}-S_{\infty}^{2}=\lim_{t\to\infty}G(t)=\lim_{t\to\infty}\int_{0}^{t}\frac{t}{y^{2}+t^{2}}\,\sin(y^{2}+t^{2})\,dy\,.
$$
We split
$$
\int_{0}^{t}\frac{t}{y^{2}+t^{2}}\,\sin(y^{2}+t^{2})\,dy=\int_{0}^{1}\frac{t}{y^{2}+t^{2}}\,\sin(y^{2}+t^{2})\,dy+\int_{1}^{t}\frac{t}{y^{2}+t^{2}}\,\sin(y^{2}+t^{2})\,dy
$$
and we observe that
$$
\left|\int_{0}^{1}\frac{t}{y^{2}+t^{2}}\,\sin(y^{2}+t^{2})\,dy\right|\le\frac{1}{t}
$$
whereas
\begin{equation*}
\begin{split}
\int_{1}^{t}\frac{t}{y^{2}+t^{2}}\,\sin(y^{2}+t^{2})\,dy&=\left[-\frac{t\cos(y^{2}+t^{2})}{2y(y^{2}+t^{2})}\right]_{y=1}^{y=t}\\
&\qquad-\frac{1}{2t}\int_{1}^{t}\frac{3t^{2}y^{2}+t^{4}}{y^{2}(y^{4}+2t^{2}y^{2}+t^{4})}\,\cos(y^{2}+t^{2})\,dy
\end{split}
\end{equation*}
and one can easily see that each term tends to zero as $t\to\infty$. Hence $(*)$ is proved.
Let us start with a warm-up exercise. Introduce the functions
$$g_{\pm}(x)=\operatorname{Ai}(x)\pm i\operatorname{Bi}(x).$$
Computing the Wronskian of these two solutions of the Airy equation, one can check that
$$\frac{1}{\operatorname{Ai}^2(x)+\operatorname{Bi}^2(x)}=\frac{\pi}{2i}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]$$
This gives the integral $\mathcal{K}(0)$ as
$$\mathcal{K}(0)=\pi\left[\arg g_+(\infty)-\arg g_+(0)\right]=\pi\left[\pi-\frac{\pi}{3}\right]=\frac{\pi^2}{6}.$$
To compute the integral $\mathcal{K}(3n)$, we will need to develop a more sophisticated approach. First note that (see here)
$$g_{\pm}(x)=-2e^{\mp 2\pi i/3}\operatorname{Ai}\left(e^{\mp2\pi i/3}x\right).$$
Therefore
\begin{align}\mathcal{K}(3n)&=\frac{\pi}{2i}\int_0^{\infty}x^{3n}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]dx=\\
&=\frac{\pi}{2i}\lim_{R\rightarrow\infty}\int_{S_R}z^{3n}\frac{\operatorname{Ai}'(z)}{\operatorname{Ai}(z)}dz,
\end{align}
where the contour $S_R$ in the complex $z$-plane
is composed of two segments: one going from $Re^{2\pi i/3}$ to $ 0$ and another one going from $0$ to $ Re^{-2\pi i/3}$.
It is a well-known fact that the Airy function $\operatorname{Ai}(z)$ has zeros (i.e. our integrand has poles) on the negative real axis only. Therefore by residue theorem our integral is equal to
$$\mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{R\rightarrow \infty}\int_{C_R}z^{3n}\left[\ln\operatorname{Ai}(z)\right]'dz,\tag{1}$$
where $C_R$ is the arc of the circle of radius $R$ centered at the origin going counterclockwise from $Re^{-2\pi i/3}$ to $Re^{2\pi i/3}$.
The limit (1), on the other hand, can be computed using the asymptotics of the Airy function as $z\rightarrow\infty$ for $|\arg z|<\pi$:
\begin{align}
\ln\operatorname{Ai}(z)\sim -\frac23 z^{3/2}-\ln2\sqrt{\pi}-\frac14\ln z+
\ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}\left(\frac43 z^{3/2}\right)^{-k}.\tag{2}
\end{align}
Note that if we introduce instead of $z$ the variable $s=\frac43z^{3/2}$, then the integration will be done over the circle of radius $\Lambda=\frac43 R^{3/2}$, i.e. a closed contour in the complex $s$-plane. The corresponding integral can therefore be computed by residues by picking the necessary term in the large $s$ expansion of $\ln \operatorname{Ai}(z)$.
More precisely, we have the following formula:
\begin{align}
\mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{\Lambda\rightarrow \infty}\oint_{|s|=\Lambda}
\left(\frac{3s}{4}\right)^{2n}d\left[-\frac16\ln s+\ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}s^{-k}\right]\tag{3}
\end{align}
To compute the residue, it suffices to expand the logarithm-sum up to order $2n$ in $s^{-1}$. Note that the Pochhammer symbol coefficients are in fact some rational numbers.
In the simplest case $n=0$, the residue is determined by the term $-\frac16\ln s$ and we readily reproduce the previous result
$$\mathcal{K}(0)=-\frac{\pi}{2i}\cdot 2\pi i\cdot\left(-\frac16\right)=\frac{\pi^2}{6}.$$
The general formula for arbitrary $n$ would look a bit complicated (but straightforward to obtain) due to the need to expand the logarithm of the sum.
Example. The calculation of the corresponding values $M(n)=\mathcal{K}(3n)$ in Mathematica can be done using the command
\begin{align}\mathtt{\text{ M[n_] := -Pi^2 SeriesCoefficient[
Series[(3 s/4)^(2 n) D[-Log[s]/6 +}} \\ \mathtt{\text{ Log[Sum[(-1)^k
Pochhammer[1/6, k] Pochhammer[5/6, k]/(k! s^k), }} \\
\mathtt{\text{{k, 0, 2 n}]], s], {s, Infinity, 1}], 1]}} \end{align}
This yields, for instance,
$$M(0)=\frac{\pi^2}{6},\quad M(1)=\frac{5\pi^2}{32},\quad M(2)=\frac{565\pi^2}{512},$$ $$\ldots, M(10)=\frac{2\,660\,774\,144\,147\,177\,521\,025\,228\,125\,\pi^2}{2\,199\,023\,255\,552},\ldots$$
and so on.
Added: The large $s$ expansion (2) can also be found directly by using Airy equation. Moreover, by transforming it into an equation for $\ln \operatorname{Ai}(z)$, one can avoid reexpanding the logarithm of the sum. The price to pay will be that the expansion coefficients will be determined by a nonlinear recurrence relation instead of explicit formulas.
Best Answer
Let $I$ denote the integral in question. With the change of variable $v = \frac{\pi x^2}{2}$, we have
$$ I = \frac{1}{\pi} \int_{0}^{\infty} \left\{ (1 - 2 C(x) )^{2} + (1 - 2S(x))^{2} \right\}^{2} \, dv $$
where $x = \sqrt{2v / \pi}$ is understood as a function of $v$. By noting that
$$ 1-2 S(x) = \sqrt{\frac{2}{\pi}} \int_{v}^{\infty} \frac{\sin u}{\sqrt{u}} \, du \quad \text{and} \quad 1-2 C(x) = \sqrt{\frac{2}{\pi}} \int_{v}^{\infty} \frac{\cos u}{\sqrt{u}} \, du, $$
we can write $I$ as
$$ I = \frac{4}{\pi^{3}} \int_{0}^{\infty} \left| A(v) \right|^{4} \, dv \tag{1} $$
where $A(v)$ denotes the function defined by
$$ A(v) = \int_{v}^{\infty} \frac{e^{iu}}{\sqrt{u}} \, du. $$
Now we want to simplify $\left| A(v) \right|^2$. To this end, we note that for $\Re u > 0$,
$$ \frac{1}{\sqrt{u}} = \frac{1}{\Gamma\left(\frac{1}{2}\right)} \frac{\Gamma\left(\frac{1}{2}\right)}{u^{1/2}} = \frac{1}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-ux}}{\sqrt{x}} \, dx = \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} e^{-ux^{2}} \, dx \tag{2} $$
Using this identity,
\begin{align*} A(v) &= \frac{2}{\sqrt{\pi}} \int_{v}^{\infty} e^{iu} \int_{0}^{\infty} e^{-u x^2} \, dx du = \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \int_{v}^{\infty} e^{-(x^2-i)u} \, du dx \\ &= \frac{2 e^{iv}}{\sqrt{\pi}} \int_{0}^{\infty} e^{-v x^2} \int_{0}^{\infty} e^{-(x^2-i)u} \, du dx = \frac{2 e^{iv}}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-v x^2}}{x^2-i} \, dx. \end{align*}
Thus by the polar coordinate change $(x, y) \mapsto (r, \theta)$ followed by the substitutions $r^2 = s$ and $\tan \theta = t$, we obtain
\begin{align*} \left| A(v) \right|^2 &= A(v) \overline{A(v)} = \frac{4}{\pi} \int_{0}^{\infty} \int_{0}^{\infty} \frac{e^{-v (x^2+y^2)}}{(x^2-i)(y^2 + i)} \, dxdy \\ &= \frac{4}{\pi} \int_{0}^{\infty} \int_{0}^{\frac{\pi}{2}} \frac{r e^{-v r^2}}{(r^2 \cos^{2}\theta-i)(r^2 \sin^{2}\theta + i)} \, d\theta dr \\ &= \frac{2}{\pi} \int_{0}^{\infty} \int_{0}^{\frac{\pi}{2}} \frac{e^{-v s}}{(s \cos^{2}\theta-i)(s \sin^{2}\theta + i)} \, d\theta ds \\ &= \frac{2}{\pi} \int_{0}^{\infty} \frac{e^{-v s}}{s} \int_{0}^{\frac{\pi}{2}} \left( \frac{1}{s \cos^{2}\theta-i} + \frac{1}{s \sin^{2}\theta + i} \right) \, d\theta ds \\ &= \frac{2}{\pi} \int_{0}^{\infty} \frac{e^{-v s}}{s} \int_{0}^{\infty} \left( \frac{1}{s -i(t^2 + 1)} + \frac{1}{s t^2 + i (t^2 + 1)} \right) \, dt ds. \end{align*}
Evaluation of the inner integral is easy, and we obtain
\begin{align*} \left| A(v) \right|^2 &= 2 \int_{0}^{\infty} \frac{e^{-v s}}{s} \Re \left( \frac{i}{\sqrt{1 + is}} \right) \, ds. \end{align*}
Applying $(2)$ again, we find that
\begin{align*} \left| A(v) \right|^2 &= 2 \int_{0}^{\infty} \frac{e^{-v s}}{s} \Re \left( \frac{i}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-(1+is)u}}{\sqrt{u}} \, du \right) \, ds \\ &= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-v s}}{s} \int_{0}^{\infty} \frac{e^{-u} \sin (su)}{\sqrt{u}} \, du\, ds \\ &= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-u}}{\sqrt{u}} \int_{0}^{\infty} \frac{\sin (su)}{s} \, e^{-v s} \, ds\, du \\ &= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-u}}{\sqrt{u}} \arctan \left( \frac{u}{v} \right) \, du \\ &= \frac{4\sqrt{v}}{\sqrt{\pi}} \int_{0}^{\infty} e^{-vx^{2}} \arctan (x^2) \, dx \qquad (u = vx^2) \tag{3} \end{align*}
Here, we exploited the identity
$$ \int_{0}^{\infty} \frac{\sin x}{x} e^{-sx} \, dx = \arctan \left(\frac{1}{s}\right), $$
which can be proved by differentiating both sides with respect to $s$.
Plugging $(3)$ to $(1)$ and applying the polar coordinate change, $I$ reduces to
\begin{align*} I &= \frac{64}{\pi^{4}} \int_{0}^{\infty} \int_{0}^{\infty} \int_{0}^{\infty} v e^{-v(x^{2}+y^{2})} \arctan (x^2) \arctan (x^2) \, dx dy dv \\ &= \frac{64}{\pi^{4}} \int_{0}^{\infty} \int_{0}^{\infty} \frac{\arctan (x^2) \arctan (x^2)}{(x^2 + y^2)^2} \, dx dy \\ &= \frac{64}{\pi^{4}} \int_{0}^{\frac{\pi}{2}} \int_{0}^{\infty} \frac{\arctan (r^2 \cos^2 \theta) \arctan (r^2 \sin^2 \theta)}{r^3} \, dr d\theta \\ &= \frac{32}{\pi^{4}} \int_{0}^{\frac{\pi}{2}} \int_{0}^{\infty} \frac{\arctan (s \cos^2 \theta) \arctan (s \sin^2 \theta)}{s^2} \, ds d\theta. \qquad (s = r^2) \tag{4} \end{align*}
Now let us denote
$$ J(u, v) = \int_{0}^{\infty} \frac{\arctan (us) \arctan (vs)}{s^2} \, ds. $$
Then a simple calculation shows that
$$ \frac{\partial^{2} J}{\partial u \partial v} J(u, v) = \int_{0}^{\infty} \frac{ds}{(u^2 s^2 + 1)(v^2 s^2 + 1)} = \frac{\pi}{2(u+v)}. $$
Indeed, both the contour integration method or the partial fraction decomposition method work here. Integrating, we have
$$ J(u, v) = \frac{\pi}{2} \left\{ (u+v) \log(u+v) - u \log u - v \log v \right\}. $$
Plugging this to $(4)$, it follows that
\begin{align*} I &= -\frac{64}{\pi^{3}} \int_{0}^{\frac{\pi}{2}} \sin^2 \theta \log \sin \theta \, d\theta = -\frac{16}{\pi^{3}} \frac{\partial \beta}{\partial z}\left( \frac{3}{2}, \frac{1}{2} \right) \end{align*}
where $\beta(z, w)$ is the beta function, satisfying the following beta function identity
$$ \beta (z, w) = 2 \int_{0}^{\infty} \sin^{2z-1}\theta \cos^{2w-1} \theta \, d\theta = \frac{\Gamma(z)\Gamma(w)}{\Gamma(z+w)}. $$
Therefore we have
\begin{align*} I &= \frac{16}{\pi^{3}} \frac{\Gamma\left(\frac{3}{2}\right)\Gamma\left(\frac{1}{2}\right)}{\Gamma(2)} \left\{ \psi_{0} (2) - \psi_{0} \left(\tfrac{3}{2} \right) \right\} = \frac{8}{\pi^2} \int_{0}^{1} \frac{\sqrt{x} - x}{1 - x} \, dx = \frac{8 (2 \log 2 - 1)}{\pi^2}, \end{align*}
where $\psi_0 (z) = \dfrac{\Gamma'(z)}{\Gamma(z)}$ is the digamma function, satisfying the following identity
$$ \psi_{0}(z+1) = -\gamma + \int_{0}^{1} \frac{1 - x^{z}}{1 - x} \, dx. $$