@Carl Love  You wrote  " I think that `+`(U[]) is faster than add(u, u= U)".  But the following examples with direct computation the sum of large number of terms does not confirm this conclusion.

The first 2 examples with the precomputed list:

restart;

X:=[seq(i^3, i=1..5000000)]:

t:=time():

`+`(X[]);

time()-t;

                            156250062500006250000000000

                                               9.453

restart;

X:=[seq(i^3, i=1..5000000)]:

t:=time():

add(x, x=X);

time()-t;

                             156250062500006250000000000

                                                7.797

In addition 2 examples without the precomputed list:

restart;

t:=time():

add(i^3, i=1..5000000);

time()-t;

                             156250062500006250000000000

                                                1.953

restart;

t:=time():

`+`(seq(i^3, i=1..5000000));

time()-t;

                             156250062500006250000000000

                                                18.032

@Alejandro Jakubi   In the help for  add  and  mul  of Maple 2015 is written:  

"The add(x) calling sequence is equivalent to add(i, i in x). Likewise, mul(x) is equivalent to mul(i, i in x)".

@Markiyan Hirnyk   Remove the extra spaces between  Matrix  and the left parenthesis from original code.

@Markiyan Hirnyk  Thank you for your helpful comment. It's interesting that  diff  command works elementwise for lists. For example, for  int  command this is not true.

@Harry Garst  Give your original pointplot, and the community will help you animate it.

@descartes0000 

expand((x+28651/50279)*(x^2+x*21628/50279+50279/28651));

evalf[20](%);

@ghoust  Just plot  f(x)  and the region under it. The area of this region is your integral:

plot(f(x), x=0..1, color=green, filled=true);

@Markiyan Hirnyk  You are wrong. The output of  convert(x/abs(x), piecewise);  is incorrect because

eval(x/abs(x), x = 0);

     Error, numeric exception: division by zero

@Markiyan Hirnyk 

Obviously OP made a mistake when calculating the derivative of  abs(x)  with Maple. Instead of his line should be

In fact, OP has calculated the second derivative of  abs(x) rather than the first derivative. I showed how to make the check what he meant.

@Markiyan Hirnyk  And if instead of  6  would be  136?

@Preben Alsholm Thank you for the deficiencies noted and useful comments.

It is likely that in this case more naturally not write a separate procedure, and simply apply the algorithm to the specific example:

restart;

Eq:=diff(y(x),x)=(y(x)-x)^2;  y(0):=0;  h:=0.1;

y[0]:=y(0);

for k to 10 do

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

od:

[seq([i*h, y[i]], i=0..10)];

@Markiyan Hirnyk  I think that for OP is required not to use ready-made solution, and write own code to the Euler method. Of course it is difficult for a beginner.

See the solution to your previous question (in the  previous thread).

@sami131   Indeed the curves intersect at about  z=11.7236908

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

Sol(11.7236908);

The default  rkf45  for initial value problems (IVP) is a Runge-Kutta Fehlberg method that produces a fifth order accurate solution.

 Hi  Earl Saltzman !

As a basis for building I took the tetrahedron with vertices [0,0,0], [2,0,0], [1,sqrt(3),0], [1,sqrt(3)/3,2*sqrt(6)/3] .  Then I find the 3 lines of intersection of the spheres below the plane xOy. The command  eliminate allows me to find the projections of these lines on the plane xOy. Therefore, the projection of the bottom part of the sphere on xOy is the region in the plane xOy bounded by the bottom line y1, and the top by the piecewise line of the lines y2 and y3. The remaining three parts  (Side1, Side2 and Side3)  are easy to get through the rotations of the Bottom .

The above code can be used to animate the body. It is only necessary to reduce  numpoints from 10000 to 2500

K := plots[display](Side1, Side2, Side3, Bottom, axes = none):

plots[animate](plottools[rotate], [K, phi, [[1, (1/3)*sqrt(3), 0], [1, (1/3)*sqrt(3), 1]]], phi = 0 .. 2*Pi);