Prove: A radial harmonic function $f$ on $\mathbb{R}^N \setminus \{0\}$ is of the form $\frac{b}{|x|^{N-2}} + c$ for $b,c \in \mathbb{R}$.
My try: Label $g_i = (0,..,0,x_i,0,..,0)$. From the maximum principle of harmonic functions and from the fact that $f$ is radial, we know that for all $i$, $f(g_i)$ either monotonically increases or decreases.($f$ can also be constant, but that's a simple case). Assume WLOG that it's monotonically increasing. Therefore: $(f(g_i))' > 0 \rightarrow f'(g_i)(0,...,0,1,0,...,0) > 0$ from the chain rule. In other words, every entry in the gradient of $f$ is positive for all $x \in \mathbb{R}^n \setminus \{0\}$. In addition, I know that form symmetry, the fact that the Laplacian is equal to $0$ means that all the seconds derivatives are $0$.
I'm not sure how to connect these things together and continue from here.
The gradient of a harmonic function is a vector field of zero divergence. Therefore (by the divergence theorem), the flux of the gradient across the sphere $S_r = \{x:|x|=r\}$ is the same for every $r$. For a radial function $f(x) = h(|x|)$, the gradient is $$\nabla f(x) = h'(|x|)\frac{x}{|x|}$$ and therefore the flux across $S_r$ is $$ \omega_{N-1} r^{N-1}h'(r) $$ where $\omega_{N-1}$ is the surface area of the unit sphere. This being constant implies that $h'(r) = cr^{1-N}$, hence $$ h(r) = A + B\cdot \begin{cases} \log r ,\quad & N=2 \\ r^{2-N},\quad & N>2 \end{cases} $$