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

@tomleslie 

m:=< <"foobar"  | 77>,

     <"faabar"  | 81>,

     <"foobaa"  | 22>,

     <"faabaa"  | 8>,

     <"faabian" | 88>,

     <"foobar"  | 27>,

     <"fiijii"  | 52>>:

keys:=["foo", "faa"]:

< seq( < j | add( `if`(m[i,1][1..3]=j,  m[i,2] ,0), i=1..op(m)[1])>, j=keys) >;

@Markiyan Hirnyk  I already wrote that the condition  (y*z+x)*y = c  is equivalent  (y/sqrt(c)*z+x/sqrt(c))*y/sqrt(c)=1 .  The objective function x+2*y*sqrt(z^2+1)  in new variables  x1=x/sqrt(c)  and  y1=y/sqrt(c)  will be  sqrt(c)*(x1+2*y1*sqrt(z^2+1))

@Jahani_21   

int(sqrt(x), x=0..4);  # Absolutely the exact result as the fraction

evalf(int(sqrt(x), x=0..4));  # The exact result as the float with 10 digits

@Carl Love  I had no idea that the solution can be so fast. Vote up!