Calculating the probability to find two bosons at the same place I want to calculate the probability to find two bosons at the same place.
Let the bosonic wave function be 
$$\Psi_B(x_1,x_2)=\frac{1}{\sqrt{2}}\left(\Psi_1(x_1)\Psi_2(x_2)+\Psi_2(x_1)\Psi_1(x_2)\right).$$
Then I have to calculate the expectation value of $\delta(x_1-x_2):$
$$\begin{align}\langle\delta(x_1-x_2)\rangle_{\Psi_B}&=\int_{\mathbb{R}^3}\int_{\mathbb{R}^3}\overline{\Psi_B(x_1,x_2)}\delta(x_1-x_2)\Psi_B(x_1,x_2)d^3x_2d^3x_1\\
&=\frac{1}{2}\int_{\mathbb{R}^6}\overline{2\Psi_1(x_1)\Psi_2(x)}\cdot2\Psi_1(x_1)\Psi_2(x)d^6x\\
&=2||\Psi_1\Psi_2||^2=2>1.\end{align}$$
This result is obviously wrong. Why is that, though?
 A: 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}
A: To calulate the the probability to find two bosons at the same place to need to evaluate 
$$\begin{align}\langle\delta(x_1-x_2)\rangle_{\Psi_B}&=\int_{\mathbb{R}^3}|{\Psi_B(x,x)}|^2d^3x\\
&=\int_{\mathbb{R}^3}|\Psi_1(x)|^2|\Psi_2(x)|^2d^3x.\end{align}$$
This is clearly smaller than 1 since you are multiplying together two normalized probability distributes and integrating over them.
