You are wrong. $x_n=O(\sqrt n$.

Let $y_n=\frac{(n-1)!!}{n!!}.$ Consider \[I_n=\int_0^{\frac{\pi}{2}}\sin^n xdx=-\cos x \sin^{n-1} x|_0^{\frac{\pi}{2}}+(n-1)(I_{n-2}-I_n)\] Therefore $I_n=\frac{n-1}{n}I_{n-2}$. Obviosli $I_0=\frac{\pi}{2}, I_1=1$. Therefore $I_n=y_n$ if $n$ - odd and $I_n=y_n\frac{\pi}{2}$ if $n$ is even. $x_n=\frac{(4n+1)!!(2n-2)!!}{(4n)!!(2n-1)!!}=\frac{(4n+1)y_{4n}}{(2n-1)y_{2n-2}}.$ $y_ny_{n-1}=\frac{1}{n}$. Consider $I_n$ for big n, then \[I_n=\int_0^{\frac{\pi}{2}}\cos^n xdx\] and $1-\frac{x^2}{2}<\cos x<1-\frac{x^2}{2}+\frac{x^4}{24}.$ Let \[I_{n0}=\int_0^{n^{-1/3}}\cos x^ndx, I_{n1}=\int_{n^{-1/3}}^{\frac{\pi}{2}}\cos^n xdx <\pi exp(-n^{1/3}/2)).\] \[{I_{n0}=\frac{1}{\sqrt n}\int_0^{n^{1/6}}exp(n\ln \cos (\frac{y}{2\sqrt n}}))dy.\] It gives $I_n=\sqrt{\frac{\pi}{2n}}(1-\frac{1}{36n}+O(\frac{1}{n^2})).$ Therefore \[x_n=\frac{4n+1}{2n-1}\sqrt{\frac{2n-2}{4n}}(1-\frac{1}{36*4n}+\frac{1}{36*2n}+O(\frac{1}{n^2}))=\sqrt 2 (1+\frac{19}{72n}+O(\frac{1}{n^2})).\]

Actually the above computation is not correct. In fact we get \[I_n=\sqrt{\frac{\pi}{2n}} \left(1-\frac{1}{4n}+\mathcal{O}\left(\frac{1}{n^2}\right)\right)\]and hence \[x_n=\sqrt{2}+\frac{5\sqrt{2}}{16n}+\mathcal{O}\left(\frac{1}{n^2}\right).\]This latter asymptotic can in fact be proved much easier: Note that \begin{align*} \log x_n&=\log\left(1+\frac{1}{2n}\right)+\log\left(1+\frac{1}{2n+2}\right)+\dots+\log\left(1+\frac{1}{4n}\right)\\ &=\left(\frac{1}{2n}+\frac{1}{2n+2}+\dots+\frac{1}{4n}\right)-\frac{1}{2}\left(\frac{1}{(2n)^2}+\frac{1}{(2n+2)^2}+\dots+\frac{1}{(4n)^2}\right)+\mathcal{O}\left(\frac{1}{n^2}\right)\\ &=\frac{1}{2n}+\frac{H_{2n}}{2}-\frac{H_n}{2}-\frac{1}{8}\left(S_{n}-S_{2n}\right)+\mathcal{O}\left(\frac{1}{n^2}\right) \end{align*}where \[H_n=1+\frac{1}{2}+\dots+\frac{1}{n}=\log n+\gamma+\frac{1}{2n}+\mathcal{O}\left(\frac{1}{n^2}\right)\]is the usual harmonic series and \[S_n=\frac{1}{n^2}+\frac{1}{(n+1)^2}+\dots=\frac{1}{n}+\mathcal{O}\left(\frac{1}{n^2}\right).\]Hence \[\log x_n=\frac{\log 2}{2}+\frac{3}{8n}-\frac{1}{16n}+\mathcal{O}\left(\frac{1}{n^2}\right)=\frac{\log 2}{2}+\frac{5}{16n}+\mathcal{O}\left(\frac{1}{n^2}\right)\]so that \[x_n=\sqrt{2} \cdot \exp\left(\frac{5}{16n}+\mathcal{O}\left(\frac{1}{n^2}\right)\right)=\sqrt{2} \cdot \left(1+\frac{5}{16n}+\mathcal{O}\left(\frac{1}{n^2}\right)\right)\]as desired.

However note that neither my nor Rust's argument actually solves the problem since we really care about an exact inequality and not an asymptotic one (we would have to make that $\mathcal{O}$ eror explicit and check cases of small $n$ by hand to make use of these estimates). Instead we consider \[y_n=\frac{(2n+2)(2n+4)\dots (4n)}{(2n+1)\dots (4n-1)}.\]Note that $x_ny_n=\frac{4n+1}{2n}=2+\frac{1}{2n}$. On the other hand, writing \[x_n=\left(1+\frac{1}{2n}\right)\left(1+\frac{1}{2n+2}\right)\dots \left(1+\frac{1}{4n}\right)\]and \[y_n=\left(1+\frac{1}{2n+1}\right)\left(1+\frac{1}{2n+3}\right)\dots \left(1+\frac{1}{4n-1}\right),\]it's easy to see that \[\left(1+\frac{1}{4n}\right)y_n<x_n<y_n \cdot \left(1+\frac{1}{2n}\right).\]Multiplying both sides by $x_n$ we find that \[\left(1+\frac{1}{4n}\right)\left(2+\frac{1}{2n}\right)<x_n^2<\left(2+\frac{1}{2n}\right)\left(1+\frac{1}{2n}\right)\]and hence \[\frac{1}{n}<x_n^2-2<\frac{3}{2n}+\frac{1}{4n^2}<\frac{2}{n}.\]In particular $x_n>\sqrt{2}$ so that \[x_n-\sqrt{2}=\frac{x_n^2-2}{x_n+\sqrt{2}}<\frac{x_n^2-2}{2\sqrt{2}}<\frac{1}{\sqrt{2}n}<\frac{2}{n}.\]The other direction is slightly more tricky as we require an upper bound for $x_n$. To conclude \[x_n-\sqrt{2}=\frac{x_n^2-2}{x_n+\sqrt{2}}>\frac{1}{(x_n+\sqrt{2})n}>\frac{1}{4n}\]we only need the weak upper bound $x_n+\sqrt{2}<4$ which follows directly from $\sqrt{2}<2$ and $x_n<2$ which is true because $x_n^2-2<\frac{2}{n} \le 2$ implies $x_n^2<4$.