sumsqint

Example:

sumsqint(1) % returns [0,1]
sumsqint(65) % returns [1,8;4,7]

sumsqint can work with large numbers if they are factored, e.g. solve x^2+y^2=(25e6)^2:
sumsqint(repmat(factor(25e6),1,2)) % returns 9x2 array

sumsqint works with variable precision integers if you have John D'Errico's Variable Precision Integer toolbox (20 July 2009 update) - see links below.

sumsqint(vpi(5)^100) % returns 51x2 vpi array
sumsqint(vpi(65)^100) % returns 5101x2 vpi array

The algorithm uses Gaussian integer factorization. See the help for more details.