If  [x,y]  is a solution then it's obviously  abs(x), abs(y)<=floor(sqrt(n))

Firstly we look for solutions  0<=x<=y , and then, using the symmetry of the equation, reproduce them:

quadsum:=proc(n::nonnegint)

local k, m, x, y, Sol;

k:=0; m:=floor(sqrt(n));

for x from 0 to m do

for y from x to m do

if x^2+y^2=n then k:=k+1; Sol[k]:=[x,y],[y,x],[-x,y],[y,-x],[-x,-y],[-y,-x],[x,-y],[-y,x]  fi;

od; od;

Sol:=convert(Sol, set);

if op(Sol)::symbol then return {} else Sol fi;

end proc:

Example:

quadsum(5000); nops(%);

Verification:

isolve(x^2+y^2=5000); nops([%]);

See  Student[MultivariateCalculus][Jacobian]

f := x->piecewise(0 < x and x < Pi, 0, Pi < x and x < 2*Pi, Pi):

a := 0: b := 2*Pi: p := b-a:

fp := f(x-floor((x-a)/p)*p):

plot([fp, seq([Pi*k, t, t = 0 .. Pi], k = -5 .. 6)], x = -6*Pi .. 6*Pi, color = [red, `$`(black, 12)], thickness = [3, `$`(1, 12)], linestyle = [1, `$`(2, 12)], discont = true, scaling = constrained);

Unfortunately Maple does not solve your equation in the general case. But it is easy to write a procedure that solves the equation   x^2+y^2 = z^2+t^2 = n   for any specified  n .  Here is the such procedure:

Sol := proc(n)

local S, L, k, i, j;

S := [isolve(x^2+y^2 = n)];

k := 0;

for i in S do

for j in S do

k := k+1; L[k] := [op(i), z = rhs(j[1]), t = rhs(j[2])]

end do end do;

L := convert(L, list);

if op(L)::symbol then return [] else L end if;

end proc:

Example of use:

S := Sol(1000):  nops(S);  seq(S[25*i], i = 1 .. 10);  # Total 256 solutions. Displayed  10  ones

restart;

combinat[permute]([0$4,1$4], 4):

L:=%[2..-1]:

P:=combinat[permute]([0$3,1$3], 3):

k:=0:

for l in L do

for p in P do

k:=k+1: M[k]:=[l,l*p[1],l*p[2],l*p[3]]:

od: od:

M:=convert(M, list):

S:={seq(seq([op(M[i,2..k-1]),M[i,1],op(M[i,k..4])], k=2..5), i=1..nops(M))}:  # Each matrix is defined as the list of lists

nops(S);

seq(Matrix(S[20*i]), i=1..10);  # As example 10 matrices from 225 ones

Edited. The code is optimized (removed unnecessary permutations)

Maybe you just want to bring your equation of order 2 to the canonical form? This is an imaginary ellipse:

eq:=25*(y1-3)^2+200+100*(x1-1)^2=0:

 (eq-200)/200;

solve(simplify(25*(y1-3)^2+200+100*(x1-1)^2=0, {(y1-3)^2+(x1-1)^2=a}), a);

 or

isolate(simplify(25*(y1-3)^2+200+100*(x1-1)^2=0, {(y1-3)^2+(x1-1)^2=a}), a);

eq:=diff(y(x), x)+4*y(x)^3-3*y(x)=0:

DEtools[DEplot](eq, y(x), x=-1 .. 0.1, y=-1 .. 0.1);

1) For the numerical solution you must specify numerical values for all parameters  (C_d, rho, r, m, g )

2) I do not understand the notation  d2v_x . If this is the second derivative, then it should be written as  diff(v_x, t,t)  

Using this procedure, you can find the values of the unknown function at certain points and build its plot:

dy4:=diff(y(x),x):

eqn4:=dy4=sin(x*y(x)):

ic1:=y(0)=1:

ans3:=dsolve({eqn4,ic1}, y(x), numeric);

ans3(1);

plots[odeplot](ans3, [x,y(x)], x=0..10, thickness=2, numpoints=1000);

The example - all on the same plot:

A := plot3d(x^2-y^2, x = -1 .. 1, y = -1 .. 1, shading = zhue):

B := plot3d([1, y, z], y = -2 .. -3/2, z = -1 .. 1, shading = zhue):

plots[display](A, B, axes = frame, orientation = [40, 75], lightmodel=light4);

map(simplify@sum, op(sum(a*u[k]+b*v[k], k=1..n)));

combine(%);

restart;

f:=(x,y)->x*(x+y)*exp(y-x);

extrema(f(x,y), {}, {x,y}, `s`);

Points:=`s`;

Student[MultivariateCalculus][SecondDerivativeTest](f(x,y), [x,y]=[0,0]);

Student[MultivariateCalculus][SecondDerivativeTest](f(x,y), [x,y]=[1/2,-3/2]);

You have forgotten  expand  f  in the Taylor series in the neighborhood  x=4. Should be

f := x^(6*ln(x));

Digits:=15;

taylor(f, x=4, 7);

T2 := convert(%, polynom);

f_value := evalf(subs(x = 5, T2));

Let   F(x,y,t,s) = K(x,y,t,s)*h_m(x)*h_n(y)*h_p(t)*h_q(s) . 

int(F(x,y,t,s), [x=0..1, y=0..1, t=0..1, s=0..1]);