$$V(x) = \frac{|y|^2-|x|^2}{|y-x|^N}$$
for fixed $y\in\mathbb{R}^n$, and $x\in\mathbb{R}^n-\{y\}$
I need to prove that this this is harmonic. That is, the sum of its second partial derivatives are $0$. Let's calculate them.
$$\left(\frac{f}{g}\right)' = \frac{fg'-gf'}{g^2}$$
so:
$$\frac{\partial V(x)}{\partial x_j} = \frac{(|y|^2-|x|^2)N|y-x|^{N-1}\frac{(-x_j)}{|y-x|} - |y-x|^N(-2|x|\frac{x_j}{|x|})}{|y-x|^{2N}} = \\\frac{-Nx_j(|y|^2-|x|^2)|y-x|^{N-2}+2x_j|y-x|^N}{|y-x|^{2N}}$$
so
$$\frac{\partial^2 V(x)}{\partial x_j^2} = \frac{\left(-Nx_j(|y|^2-|x|^2)|y-x|^{N-2}+2x_j|y-x|^N\right)\left(2N|y-x|^{2N-1}\frac{(-x_j)}{|y-x|}\right)-|y-x|^{2N}\left(\cdots\right)}{|y-x|^{4N^2}}$$
where $\left(\cdots\right)$ is the derivative of $\left(-Nx_j(|y|^2-|x|^2)|y-x|^{N-2}+2x_j|y-x|^N\right)$. It has a product of three things, and another term. It's huge. Multiplied by the other things, it's even bigger.
I think this is not the right way to solve this exercise. Is there a better way? If not, how can I do these giant calculations by hand? It's too already in the computer using LaTeX.
Use the product rule for Laplacian: $$\Delta(fg) = f\Delta g+g\Delta f+2\nabla f\cdot\nabla g\tag0$$
Here $f(x)=|x|^2-|y|^2$ and $g(x)=|x-y|^{-N}$. (I flipped the sign of $V$ to have fewer sign changes around.) So,
$$\nabla f = 2x\tag1$$ $$\nabla g = -N(x-y)|x-y|^{-N-2} \tag2$$ $$\Delta f = \operatorname{div}\nabla f = 2N \tag3$$ $$\Delta g = \operatorname{div}\nabla g = -N^2 |x-y|^{-N-2} + N(N+2) (x-y)\cdot (x-y) |x-y|^{-N-4} = 2N|x-y|^{-N-2} \tag4$$ where the last computation used the product rule $\operatorname{div}(\phi F)=\phi\operatorname{div}( F) + \nabla \phi\cdot F$.
Stick (1)-(4) into (0): $$ \Delta(fg) = 2N(|x|^2-|y|^2) |x-y|^{-N-2} + 2N|x-y|^{-N} -4N x\cdot(x-y)|x-y|^{-N-2} \tag5 $$ The rest is algebra. Let's put $z=x-y$, so $x=y+z$, and expand $|x|^2 = |y|^2 + 2y\cdot z+|z|^2$. Then (5) becomes $$ \Delta(fg) = 2N(2y\cdot z + |z|^2) |z|^{-N-2} + 2N|z|^{-N} -4N (y\cdot z + |z|^2) |z|^{-N-2} \tag6 $$ which happily collapses into $0$.