Consider the function which is defined by the equation and its symmetries:

restart;

f:=x->sin(2*x)^2-cos(8*x)^2-(1/2)*cos(10*x):

is(f(Pi-x)=f(Pi+x)) and is(f(Pi/2-x)=f(Pi/2+x));

                                        true

Therefore it is sufficient to solve this equation in the range  0 .. Pi/2

RootFinding[Analytic](f(x), x, re = 0 .. (1/2)*Pi, im = -1 .. 1);

identify([%]);

You can enter a variable (the counter), first assign it the value  0 , and then at each step of the loop is increased by 1 .

Example: find all Pythagorean triangles with hypotenuse not exceeding 100 . Solution in the for loop:

restart;

k:=0: # The counter for the steps of the loop

n:=0: # The counter for solutions

for c from 2 to 100 do

for a from 1 to c-1 do

k:=k+1:

b:=sqrt(c^2-a^2):

if type(b,integer) and b>a then n:=n+1; L[n]:=[a,b,c] fi;

od:

od:

L:=convert(L, list):

k;  n;  L;

Addition:  The example of a purely illustrative nature. There are far more effective methods of generating Pythagorean triples.

restart;

 L:=[seq(seq(x^i*y^j, j=0..3-i), i=0..3)];

 add(a[k]*L[k+1], k=0..nops(L)-1);

The procedure  Euler  solves your problem:

restart;

Euler:=proc(Eq,Ivp,h,n)  # n is the number of steps

local Y, k;

Y[0]:=Ivp[2];

for k to n do

Y[k]:=Y[k-1]+h*eval(rhs(Eq),{x=Ivp[1]+h*(k-1), y(x)=Y[k-1]});

od;

[seq([Ivp[1]+i*h,Y[i]], i=0..n)];  # List of points of the approximate solution

end proc:

Example of use:

Eq:=diff(y(x),x)=(y(x)-x)^2; Ivp:=[0,0];

Euler(Eq, Ivp, 0.1, 10);

plot([%, x-tanh(x)], x=0..1, color=[red, blue]); # Graphs of approximate (red) and exact (blue) solutions

The code of the procedure was edited.

First get more awkward answer, but after a few conversions manages to simplify it.

restart;

dsolve({diff(y(x),x)=(y(x)-x)^2, y(0)=0}, y(x));

assign(%):

convert(expand(simplify(convert(y(x), trig))), tanh);

restart;

A:=(n,a,b,c,p)->((1+c*p)*a/(1-p)+(2+c*p)*b)*mul((1+c*(1-p)^k)*p, k=1..n-1)+add((a+b+k*b*(1-p))*mul((1+c*(1-p)^j)*p,j=k..n-1), k=2..n);

for n do

a[n]:=A(n,0.2,0.09,0.57,0.3);

if a[n]>=1  then Sol:=[[n-1,a[n-1]], [n,a[n]]]; break; fi;

od:

Sol;

You choose what  answer to you is more appropriate  n=6  or  n=7

The code is edited.

The simple procedure A  constructs the characteristic matrix with your properties of any size. For example displayed  A(10). Then we find the characteristic polynomial for  n = 50, and is set to 0 the coefficient of the imaginary part. We find that it is equal to  0  only for  x = 0. Thus we have a symmetrical real matrix that have only real roots:

A:=n->Matrix(n, {(1,1)=-I*x-lambda, (n,n)=-I*x-lambda, seq((i,i)=-lambda, i=2..n-1), seq((i,i-1)=-1, i=2..n), seq((i-1,i)=-1, i=2..n)}):

A(10);

B:=sort(LinearAlgebra[Determinant](A(50)), lambda);

x=solve(coeff(B,I)=0, x);

Verification:

C:=subs(x=0, B):

fsolve(C=0, lambda, complex);

We got only real roots.

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