v := <1, 0, 0>:  u := beta -> <0, cos(beta), sin(beta)>:  c := <0, 0, 0>:  a := 2:  b := 3:

plot3d(c+a*v*sec(alpha)+b*u(beta)*tan(alpha), beta = 0 .. 2*Pi, alpha = -(1/2)*Pi .. (1/2)*Pi, style = surface, color = khaki, view = [-20 .. 20, -20 .. 20, -20 .. 20], numpoints = 50000, axes = frame);

f:=x->piecewise(-2<=x and x<0, -x^2 ,0<=x and x<=2,x^2, undefined);

solve(f(2)-f(-2)=D(f)(c)*(2-(-2)), {c});

plot([f(x), 2*x, f(-1)+D(f)(-1)*(x-(-1)), f(1)+D(f)(1)*(x-1)], x=-2..2, color=[red,black,blue,blue], linestyle=[1,2,1,1], thickness=[2,1,1,1], scaling=constrained);  # Visualization

Addition.  Maybe the mean value theorem for definite integral  is meant. 

int(f(x), x=a..b) = f(c)*(b-a) -> f(c)=int(f(x), x=a..b)/(b-a)

int(f(x), x=-2..2)/(2-(-2));

                      0

I think there is no need to use Student:-Calculus1:-RiemannSum command. You can simply use  standard  add  command:

a:=0:  b:=4:  N:=100:  f:=x->sqrt(x):  h:=(b-a)/N:

add(f(a+(n-1)*h)*h, n=1..N):

Approx:=evalf(%); Exact:=evalf(int(sqrt(x), x=0..4));

For your integral much more accurate result gives the middle Riemann sum:

a:=0:  b:=4:  N:=100:  f:=x->sqrt(x):  h:=(b-a)/N:

add(f(a+h/2+(n-1)*h)*h, n=1..N):

Approx:=evalf(%); Exact:=evalf(int(sqrt(x), x=0..4));

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(%);