Question about parametric resolution of a nonlinear system.

edit

Question about parametric resolution of a nonlinear system.

Given the set of all Pythagorean quadruples a^2+b^2+c^2=d^2

a=u^2+v^2-s^2-t^2

b=2(ut+v*s)

c=2(vt-u*s)

d=u^2+v^2+s^2+t^2

and given a subset of Pythagorean quadruples that have the characteristic d+a=c^2

d+a=c^2

d=36m^2+18m+4n^2+2n+3

a=24mn+6m+6n+1=3(3N-3)/6+1

b=2(3m+n+1)(6m-2*n+1)

c=2(3m+n+1)

to learn more about this type of Pythagorean quadruple look here

https://drive.google.com/file/d/1AgSi...

the demonstration here

https://drive.google.com/file/d/11zU-...

I don't know how to solve the system in parametric solutions but this way doesn't work

import time
Start_Time = time.time()
var('N a b c d u v s t m n')

eq0 = N -65 == 0
eq1 = u^2+v^2-s^2-t^2 - a == 0
eq2 = 2*(u*t+v*s) - b == 0
eq3 = 2*(v*t-u*s) - c == 0
eq4 = u^2+v^2+s^2+t^2 - d == 0
eq5 = d+a - c^2 == 0
eq6 = 36*m^2+18*m+4*n^2+2*n+3 - d  == 0
eq7 = 24*m*n+6*m+6*n+1 - a == 0
eq8 = 2*(3*m+n+1)*(6*m-2*n+1) - b == 0
eq9 = 2*(3*m+n+1) - c == 0
eq10 = 3*(3*N-3)/6+1 - a == 0

solutions = solve([eq0,eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10],N,a,b,c,d,u,v,s,t,m,n)
sol = solutions
Execution_Time = time.time() - Start_Time
print (Execution_Time)
print(sol)

Is there any other working method?