@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!

@Stavros   map(Conv, S);

@veg_nw  I do not know this.

@Stavros   Here is the simple procedure which converts any matrix or list of lists to this form:

Conv:=proc(A::{Matrix,listlist})

local m, n, L, i, j;

if A::Matrix then m,n:=op(A)[1],op(A)[2] else m:=nops(A); n:=nops(A[1]) fi;

L:=0;

for i to m do

for j to n do

L:=L+A[i,j]*cat(E,i,j)

od; od;

L;

end proc:

Example of work:

Conv([[0,1,0,1], [0,1,0,1], [0,0,0,0], [0,0,0,0]]);

                             E12+E14+E22+E24

@Axel Vogt   

The idea of the above code is based on the fact that if the rank of the matrix is 1, then there exists a non-zero row and every other row is proportional to it (in this case with coefficients 0 or 1).
The code can easily be generalized to the procedure for the matrix  n x n :

restart;

G:=proc(n)

local L, P, k, M, P1, l, p;

combinat[permute]([0$n,1$n], n):

L:=%[2..-1]:

P:=combinat[permute]([0$(n-1),1$(n-1)], n-1):

k:=0:

for l in L do

for p in P do

k:=k+1: M[k]:=[l,seq(l*p[i], i=1..n-1)]:

od: od:

M:=convert(M, list):

{seq(seq([op(M[i,2..k-1]),M[i,1],op(M[i,k..n])], k=2..n+1), i=1..nops(M))}:

end proc:

@CakePie  

simplify(25*(y1-3)^2+200+100*(x1-1)^2=0, {(y1-3)^2+(x1-1)^2=a});

isolate(%, a);

@Carl Love  I thought about this problem (arbitrary number nested loops), but did not know how to solve it.

To solve the problem you just can use Carl's procedure  Max :

Max:= proc(L::list)

local i, Max:= -infinity;

     for i to nops(L) do

          if L[i] > Max then

               Max:= L[i]

          end if

     end do;

     return Max

end proc:

R:=rand(1..50):

L:=[seq(R(), i=1..20)];

Max(L);

                          L := [26, 29, 29, 14, 50, 47, 8, 42, 46, 44, 9, 13, 16, 1, 6, 13, 46, 39, 28, 7]

                                                                                 50