Proving positive definiteness of matrix $a_{ij}=\frac{2x_ix_j}{x_i + x_j}$ I'm trying to prove that the matrix with entries $\left\{\frac{2x_ix_j}{x_i + x_j}\right\}_{ij}$ is positive definite for all n, where n is the number of rows/columns. 
I was able to prove it for the 2x2 case by showing the determinant is always positive. However, once I extend it to the 3x3 case I run into trouble. I found a question here
whose chosen answer gave a condition for positive definiteness of the extended matrix, and after evaluating the condition and maximizing it via software, the inequality turned out to hold indeed, but I just can't show it.
Furthermore, it would be way more complicated when I go to 4x4 and higher. I think I should somehow use induction here to show it for all n, but I think I'm missing something. Any help is appreciated.
Edit: Actually the mistake is mine, turns out there are indeed no squares in the denominator, so it turns out user141614's first answer is what I really needed. Thanks a lot! Should I just accept this answer, or should it be changed back to his first answer then I accept it?

(Update: Some fixes have been added because I solved the problem with $a_{ij}=\frac{2x_ix_j}{x_i+x_j}$ instead of $a_{ij}=\frac{2x_ix_j}{x_i^2+x_j^2}$. Thanks to Paata Ivanisvili for his comment.)
The trick is writing the expression $\frac{1}{x^2+y^2}$ as an integral.
For every nonzero vector $(u_1,\ldots,u_n)$ of reals,
$$
(u_1,\ldots,u_n) \begin{pmatrix}
a_{11} & \dots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \dots & a_{nn} \\
\end{pmatrix}
\begin{pmatrix} u_1\\ \vdots \\ u_n\end{pmatrix} =
\sum_{i=1}^n \sum_{j=1}^n a_{ij} u_i u_j =
$$
$$
=2 \sum_{i=1}^n \sum_{j=1}^n u_i u_j x_i x_j \frac1{x_i^2+x_j^2} =
2\sum_{i=1}^n \sum_{j=1}^n u_i u_j x_i x_j \int_{t=0}^\infty 
\exp \Big(-(x_i^2+x_j^2)t\Big) \mathrm{d}t =
$$
$$
=2\int_{t=0}^\infty 
\left(\sum_{i=1}^n u_i x_i \exp\big(-{x_i^2}t\big)\right)
\left(\sum_{j=1}^n u_j x_j \exp\big(-{x_j^2}t\big)\right)
\mathrm{d}t =
$$
$$
=2\int_{t=0}^\infty 
\left(\sum_{i=1}^n u_i x_i \exp\big(-{x_i^2t}\big)\right)^2
\mathrm{d}t \ge 0.
$$
If $|x_1|,\ldots,|x_n|$ are distinct and nonzero then the last integrand cannot be constant zero: for large $t$ the minimal $|x_i|$ with $u_i\ne0$ determines the magnitude order. So the integral is strictly positive.
If there are equal values among $|x_1|,\ldots,|x_n|$ then the integral may vanish. Accordingly, the corresponding rows of the matrix are equal or the negative of each other, so the matrix is only positive semidefinite.

You can find many variants of this inequality. See 
problem A.477. of the KöMaL magazine.
The entry $a_{ij}$ can be replaced by $\left(\frac{2x_ix_j}{x_i+x_j}\right)^c$ with an arbitrary fixed positive $c$; see
problem A.493. of KöMaL.
