restart;

eliminate({(x-1)^2+(y-sqrt(3))^2+z^2 = 4, (x-1)^2+(y-(1/3)*sqrt(3))^2+(z-2*sqrt(6)*(1/3))^2 = 4}, z):

y1 := solve(op(%[2]), y)[2]:

eliminate({(x-2)^2+y^2+z^2 = 4, (x-1)^2+(y-(1/3)*sqrt(3))^2+(z-2*sqrt(6)*(1/3))^2 = 4}, z):

y2 := solve(op(%[2]), y)[1]:

eliminate({x^2+y^2+z^2 = 4, (x-1)^2+(y-(1/3)*sqrt(3))^2+(z-2*sqrt(6)*(1/3))^2 = 4}, z):

y3 := solve(op(%[2]), y)[1]:

Seq := 2*sqrt(6)*(1/3)-sqrt(4-(x-1)^2-(y-(1/3)*sqrt(3))^2), x = 0 .. 2, y = y1 .. piecewise(`and`(x >= 0, x <= 1), y2, `and`(x >= 1, x <= 2), y3), style = surface, numpoints = 10000, scaling = constrained:

A := plot3d(Seq, color = red):

B := plot3d(Seq, color = blue):

C := plot3d(Seq, color = green):

Bottom := plot3d(Seq, color = yellow):

Side1 := plottools[rotate](A, 2*Pi*(1/3), [[0, 0, 0], [1, (1/3)*sqrt(3), (1/6)*sqrt(6)]]):

Side2 := plottools[rotate](B, 4*Pi*(1/3), [[0, 0, 0], [1, (1/3)*sqrt(3), (1/6)*sqrt(6)]]):

Side3 := plottools[rotate](C, (-2*Pi)*(1/3), [[2, 0, 0], [1, (1/3)*sqrt(3), (1/6)*sqrt(6)]]):

plots[display](Side1, Side2, Side3, Bottom, axes = none, orientation = [80, 80, -160]);

Maple can solve your system only numerically:

Sys:=-diff(p(z),z)=(0.855*(1+2*x(z)))/p(z), diff(x(z),z)=((6.39*10^(-3)*p(z)*(1-x(z)))/(1+2*x(z)));

Inc:=p(0)=5, x(0)=0;

Sol:=dsolve({Sys, Inc}, numeric);

plots[odeplot](Sol, [[z,p(z)], [z,x(z)]], z=0..10, color=[red,blue], thickness=2);

Addition: using the procedure  Sol , you can not only build graphics of solutions as above, but also to find solutions values at certain points, such as

Sol(5);

                           [z = 5., p(z) = 3.91451233089194739, x(z) = .120154610510706167]

1) Let  0<=a<1 .

For such  a  the limit exists and equals  0. The proof follows from well-known inequalities  sin(a)<=a  and  (a+b)/2>=sqrt(a*b)  for non-negative numbers and simple estimates

(1-cos(x*y^3))/(x^2+y^6)^(1+a) <= 2*sin(x*y^3/2)^2/(2*sqrt(x^2*y^6))^(1+a)

<=2*x^2*y^6/4/(2^(1+a)*x^(1+a)*y^(3+3*a)) = 1/2^(2+a)*x^(1-a)*y^(3-3*a)

2) If  a = 1, then the limit does not exist, as along different paths we get different results

limit(eval((1-cos(x*y^3))/(y^6+x^2)^2, x = y^3), y = 0);

limit(eval((1-cos(x*y^3))/(y^6+x^2)^2, x = y), y = 0);

The command  extrema  does not find any local or global extremes, but simply returns the maximum and minimum values of the function at the critical points found. This is illustrated by a simple example

f:=x^2-y^2:

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

This is a saddle point, there is no extremum in it.

Therefore, critical points results require a separate study, for example, using the second derivative.

We consider only the case c> 0. For simplicity, we assume  c = 1, because the general case reduces to scaling variables  x  and  y .  The condition  (y*z+x)*y = c  is equivalent  (y/sqrt(c)*z+x/sqrt(c))*y/sqrt(c)=1

F:=subs(x=solve((y*z+x)*y-1, x), x+2*y*sqrt(z^2+1));

extrema(F, {}, {y,z}, 's');  s;

Student[MultivariateCalculus][SecondDerivativeTest](F, [y,z]=[rhs(s[1,1]),rhs(s[1,2])]);

The second critical point can be studied similarly.

A := (lambda1, lambda2) -> Matrix(2, [1 - lambda1, 2, 3, 4 - lambda2]);

Example:

A(2, 4);

 Addition:  if  lambda1=lambda2  then even simplier  A := lambda -> Matrix(2, [1 - lambda, 2, 3, 4 - lambda]); 

Then

A(2)   or   A(4)

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