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

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