This problem looks hard, and may have some issue.

The overall arrangement is to an extent telescoping, since if 

a^2+b^2=c^2, then we can write 

c^2+d^2=e^2, and 

e^2+f^2=g^2, and 

g^2+h^2=j^2,

...

in each case a triple of 2 sides plus a hypotenuse.

The hypotenuse of each triple is a side of the next and so must be the odd side. This leads to the important question: how big is j?

Let v=(u-1) for its nice properties, and apply Euclid's formula iteratively, varying U,V according to context.

Then (u^2-v^2)=(2u-1) implies 2uv=(2u^2-2u), and (u^2+v^2)=(2u^2-2u+1) (a,b,c)

And 2(u^2-u+1)-1 implies 2UV=(2(u^2-u+1)^2-2(u^2-u+1))=2(u^2-u)^2 + 2(u^2-u), and (U^2+V^2)=(2(u^2-u+1)^2-2(u^2-u+1)+1)=2(u^2-u)^2 + 2(u^2-u)+1 (c,d,e)

And 2((u^2-u)^2+(u^2-u)+1)-1 implies 2UV=(2((u^2-u)^2+(u^2-u)+1)^2-2((u^2-u)^2+(u^2-u)+1))=2(u^2+(u^2-u)^2-u)^2+2(u^2+(u^2-u)^2-u) and 

(U^2+V^2)=(2((u^2-u)^2+(u^2-u)+1)^2-2((u^2-u)^2+(u^2-u)+1)+1)=2(u^2+(u^2-u)^2-u)^2+2(u^2+(u^2 - u)^2-u)+1 (e,f,g)

And 2((u-1)^2u^2+u^2+(u^2+(u^2-u)^2-u)^2-u+1)-1 (I) implies

2UV=(2((u-1)^2u^2+u^2+(u^2+(u^2-u)^2-u)^2-u+1)^2-2((u-1)^2u^2+u^2+(u^2+(u^2-u)^2-u)^2-u+1))=2(u^2+(u^2-u)^2+(u^2+(u^2-u)^2-u)^2-u)^2+2(u^2+ (u^2-u)^2+(u^2+(u^2-u)^2-u)^2-u), 

and (U^2+V^2)

=(2((u-1)^2u^2+u^2+(u^2+(u^2-u)^2-u)^2-u+1)^2-2((u-1)^2u^2+u^2+(u^2+(u^2-u)^2-u)^2-u+1)+1)

=2(u^2+(u^2-u)^2+(u^2+(u^2-u)^2-u)^2-u)^2+2(u^2+(u^2-u)^2+(u^2+(u^2-u)^2-u)^2-u)+1 (II) (g,h,j)

Setting II with (2u-1)=1499 implies that j is around 2*10^46 (before squaring) so j^2 could have as many as 93 decimal digits.

Since 5*13*17*25*29*37*41*53*61*65*73*85*89*97*101*109*113*125*137*149*157*169*173*181 or 36079677936912233935019538562910577831640625 with a mere 43 decimal digits has 2^24 or 16777216 representations as a sum of two squares, with u around 504, the maximal number of solutions for some a in the given range could be well over a million.