How to prove $$\sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{m^2+n^2}=+\infty.$$
I try to do like $$\sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{m^2+n^2}=\sum_{N=1}^\infty \sum_{n+m=N}^\infty \frac{1}{m^2+n^2}=\sum_{N=1}^\infty \sum_{m=1}^{N-1} \frac{1}{m^2+(N-m)^2}$$ $$\frac{1}{m^2+(N-m)^2}\leq \frac{2}{N^2}$$ but it doesn't work.
$$ \begin{align} \sum_{m=1}^\infty\sum_{n=1}^\infty\frac1{m^2+n^2} &\ge\sum_{m=1}^\infty\sum_{n=1}^\infty\frac1{m^2+n^2+2mn+m+n}\\ &=\sum_{m=1}^\infty\sum_{n=1}^\infty\left(\frac1{m+n}-\frac1{m+n+1}\right)\\ &=\sum_{m=1}^\infty\frac1{m+1}\\[6pt] &=\infty \end{align} $$
If the sum were finite, then we could get a contradiction as follows. Breaking it up into 4 sums depending on whether or not $m$ and $n$ are even, we have
$$\sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{m^2+n^2} = \sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{(2m)^2+(2n)^2} + \sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{(2m-1)^2+(2n)^2} $$ $$+ \sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{(2m)^2+(2n-1)^2} + \sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{(2m-1)^2+(2n-1)^2} $$ Note that each of the last three sums is greater than the first due to the denominators of each term being smaller. Thus we have $$\sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{m^2+n^2} > 4 \sum_{m=1}^\infty \sum_{n=1}^\infty \frac{1}{(2m)^2+(2n)^2} $$ But factoring out the 4 from the denominator we see the right hand side is the same as the left. Hence we get a contradiction and the sum is not finite.
You're almost there.
$$\frac{1}{m^2 + (N-m)^2} = \frac{1}{2m^2 +N^2 -2mN} = \frac{1}{2m(m-N) +N^2}\ge \frac{1}{N^2}$$
Now $$\sum_{N=1}^\infty \sum_{m=1}^{N-1}\frac{1}{m^2 + (N-m)^2} \ge \sum_{N=1}^\infty \sum_{m=1}^{N-1}\frac{1}{N^2} = \sum_{N=1}^{\infty}\frac{N-1}{N^2}$$
Can you finish from here?
$$ \begin{align} \sum_{m=1}^\infty\sum_{n=1}^m\frac1{m^2+n^2} &\ge\sum_{m=1}^\infty\sum_{n=1}^m\frac1{2m^2}\\ &=\frac12\sum_{m=1}^\infty\frac1m\\ &=+\infty. \end{align} $$
$$\iint_{(1,+\infty)^2}\frac{dx\,dy}{x^2+y^2}\geq \int_{\varepsilon}^{\pi/2-\varepsilon}\int_{\rho_0}^{+\infty}\frac{1}{\rho}\,d\rho\,d\theta = +\infty $$ hence your series is divergent by comparison with a divergent integral.
As an alternative, let
$$ r_2(n) = \left|\{(a,b)\in\mathbb{Z}^2:a^2+b^2=n\}\right| $$
and $\mathbb{1}_S(n)$ the characteristic function of the non-zero integers which are sums of two squares.
By summation by parts and by Gauss circle problem
$$ \sum_{n=1}^{N}\frac{\mathbb{1}_S(n) r_2(n)}{n}=\pi+\sum_{n=1}^{N}\frac{\pi}{N}+O(1)=\pi\log(N)+O(1).$$
$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
Note that
\begin{align} \sum_{n = 1}^{\infty}{1 \over m^{2} + n^{2}} & = \sum_{n = 1}^{\infty}{1 \over \pars{n - \ic m}\pars{n + \ic m}} = {\Psi\pars{-\ic m} - \Psi\pars{\ic m} \over -2\ic m} \\[5mm] & = {1 \over m}\,\Im\Psi\pars{\ic m}\qquad\pars{~\Psi:\ Digamma\ Function~} \\[5mm] & = {1 \over 2m^{2}} + {\pars{\pi/2}\coth\pars{\pi m} \over m} \,\,\,\stackrel{\mrm{as}\ m\ \to\ \infty}{\sim}\,\,\,{\pi \over 2}\,{1 \over m} \end{align}
such that the double series diverges because $\ds{\sum_{m = 1}^{M}{\pi \over 2}\,{1 \over m} = {\pi \over 2}\,H_{M}}$ where $\ds{H_{z}}$ is the Harmonic Number.
