How can I find the possible values that $\gcd(a+b,a^2+b^2)$ can take, if $\gcd(a,b)=1$

We have $a^2 + b^2 - a(a+b) = b^2 - ab = -b (a-b)$ and $a^2 + b^2 - b(a+b) = a^2 - ba = a (a-b)$.

So if $d$ divides both $a+b$ and $a^2+b^2$, then $d$ divides $$\gcd(a (a-b), b (a-b)) = \gcd(a, b) (a-b) = a - b.$$

So $d$ divides $a+b + a - b = 2a$ and $a+b - (a - b) = 2b$.

So $d$ divides $2\gcd(a,b)=2$.

So the possibilities for the $\gcd$ appear to be $1$ and $2$, and both clearly occur.

$$\begin{align} \gcd(a+b, a^2 + b^2) &= \gcd(a+b, a^2 + b^2 - a(a+b)) \\&= \gcd(a+b, b^2 - ab) \\&= \gcd(a+b, b^2 - ab + b(a+b)) \\&= \gcd(a+b, 2b^2) \end{align} $$

Now, $\gcd(a+b,b) = \gcd(a,b) = 1$, so we can get rid of the factors of $b$ and have

$$ \gcd(a+b, a^2 + b^2) = \gcd(a+b, 2) $$

The strategy I used was still the basic idea of the Euclidean algorithm; since I couldn't compare numeric values, I instead simplified by working to eliminate the variable $a$, starting with the largest power of $a$.

Since $\gcd(a,b)=1$, Bezout's Identity says we have an $x$ and $y$ so that $$ ax+by=1\tag{1} $$ Note that $$ \begin{align} 2a^2&=(a^2+b^2)+(a+b)(a-b)\\ 2ab&=(a+b)^2-(a^2+b^2)\\ 2b^2&=(a^2+b^2)-(a+b)(a-b)\\ \end{align}\tag{2} $$ Therefore, incorporating $(1)$ and $(2)$, $$ \begin{align} 2 &=2(ax+by)^2\\ &=2a^2x^2+4abxy+2b^2y^2\\ &=\Big((a^2+b^2)+(a+b)(a-b)\Big)x^2\\ &+2\Big((a+b)^2-(a^2+b^2)\Big)xy\\ &+\Big((a^2+b^2)-(a+b)(a-b)\Big)y^2\\ &=\color{#00A000}{(x-y)^2}\color{#C00000}{(a^2+b^2)} +\color{#00A000}{((x^2-y^2)(a-b)+2xy(a+b))}\color{#C00000}{(a+b)}\tag{3} \end{align} $$ Equation $(3)$ says that $$ \gcd(a+b,a^2+b^2)\,|\,2\tag{4} $$ Note that $\gcd(1+2,1^2+2^2)=1$ and $\gcd(1+3,1^2+3^2)=2$, so both $1$ and $2$ are possible.