If the factorization of the integer $n$ is either $n=p_1^{a_1}\cdots p_s^{a_s}$ or $n=2p_1^{a_1}\cdots p_s^{a_s}$, where all $p_i$ are primes of form $4k+1$, then $n$ has $2^{s-1}$ representations as $n=a^2+b^s$, where $\gcd(a,b)=1$ (see for example this post ).
The SAGE command sum_of_k_squares(2,n)
will give one such representation. Is there a way to find all of them?