This probability is actually $0$. This is a 2d variant of the question: "What is the probability of find the system described by $\psi(x)$ at position $x_0$?". In both cases the answer is $0$. In 1d this is because the probability is the area under a single point of height $\vert \psi(x_0)\vert^2$, and this point has no width so the area is $0$. In the 2d case of this question, the answer if $0$ because this probability is a line of finite length but no width, so the area under this line is $0$.
This does not mean that finding the expectation value of $\delta(x_1-x_2)$ cannot be done. It's just that $\langle\delta(x_1-x_2)\rangle$ does not give the probability you are looking for.
The same situation occurs in 1d. $\langle \delta (x)\rangle$ simply returns $\vert\psi(0)\vert^2$, which incidentally can be $>1$: there is nothing to prevent $\psi(x)$ or $\vert\psi(x)\vert^2$ to be $>1$ since the probability is the probability density $\vert\psi(x)\vert^2$ integrated over an interval. For small intervals $dx$ we have $\vert\psi(x)\vert^2dx<1$. An easy example of this is the normalized wavefunction $\psi(x)=\sqrt{2}\frac{e^{-2x^2}}{\pi^{1/4}}$. Then
$$
\langle \delta(0)\rangle = \int dx\, \frac{2e^{-4x^2}}{\sqrt{\pi}}\delta(0)= \frac{2}{\sqrt{\pi}}\approx 1.13.
$$
Now, in the 2d case, to find $\langle\delta(x_1-x_2)\rangle$, you obtained quite correctly the expectation value of $\delta(x_1-x_2)$ as
\begin{align}
\langle\delta(x_1-x_2)\rangle&=
\int dx_1 \int dx_2 \delta(x_1-x_2) \textstyle\frac{1}{2}
\vert \Psi_1(x_1)\Psi_2(x_2)+\Psi_2(x_1)\Psi_1(x_2)\vert^2\, ,\\
&=\int dx_1 \textstyle\frac{1}{2}\vert \Psi_1(x_1)\Psi_2(x_1)+\Psi_2(x_1)\Psi_1(x_1)\vert^2\, ,\\
\end{align}
At this point, the computation continues as
\begin{align}
&=\int dx_1 \textstyle\frac{1}{2}\vert \Psi_1(x_1)\Psi_2(x_1)+\Psi_2(x_1)\Psi_1(x_1)\vert^2\, ,\\
&=\int dx_1 \textstyle\frac{1}{2}
\vert2 \Psi_1(x_1)\Psi_2(x_1)\vert^2\, ,\\
&=2\int dx_1
\vert \Psi_1(x_1)\Psi_2(x_1)\vert^2\, ,\\
&=2\int dx_1
\vert \Psi_1(x_1)\vert^2\vert \Psi_2(x_1)\vert^2 \, , \tag{1}\\
&\ne 2 \int dx_1
\vert \Psi_1(x_1)\vert^2\int dx_1 \vert \Psi_2(x_1)\vert^2\, .
\end{align}
Of course the factor of 2 in Eq.(1) is expected since bosons are more likely to bunch together than independent particles, in other words the probability amplitude $\Psi(x,x)$ for two bosons at the same place is enhanced by $\sqrt{2}$ over the probability amplitude for distinguishable particles.
As a calculational example consider $\psi_0(x)$ and $\psi_1(x)$ to be harmonic oscillator wavefunctions:
$$
\psi_0(x)=\frac{e^{-\frac{x^2}{2}}}{\sqrt[4]{\pi }}\, ,\qquad
\psi_1(x)=\frac{\sqrt{2} e^{-\frac{x^2}{2}} x}{\sqrt[4]{\pi }}
$$
The probability density is then
$$
\vert\Psi_B(x,y)\vert^2=\frac{e^{-x^2-y^2} \left(-2 x^2 y-x \left(2 y^2-1\right)+y\right)^2}{2 \pi }
$$
Setting $x=y$ we obtain
$$
\vert\Psi_B(x,x)\vert^2=\frac{2 e^{-2 x^2} x^2 \left(1-2 x^2\right)^2}{\pi }
$$
and
\begin{align}
\int dx \vert\Psi_B(x,x)\vert^2&=\frac{7}{8\sqrt{2\pi}}\approx 0.35\, ,\\
&=2\int dx \vert\psi_0(x)\vert^2 \vert\psi_1(x)\vert^2
\end{align}
