@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

@tomleslie  y<=abs(x)  is not enough. Should be  abs(y)<=abs(x)

@Mac Dude  Of course, you're right. But there are many cases when the Maple can calculate exactly, but does not make it right. Here is an example with the angle of 27 degrees = 3*Pi/20 radians :

sin(3*Pi/20);

convert(%, radical);

On this occasion see my post

@UndergradMaths   Fastest way to input and solve a differential equation is usage  the prime marks for the derivatives in document mode:

@Carl Love  and  @Markiyan Hirnyk 

My plot is not a fake and just the extensions of the domain of the function given OP for each value of  phi . In these extensions for each of the four plots I defined the functions as a linear functions so that they are as close as possible to the original plot, specified by OP.

@Carl Love   

plots:-animate(plots:-tubeplot,  [[cos(t+phi), sin(t+phi), t], t= 0..6*Pi, radius=1/3, numpoints=200, scaling=constrained], phi= 0..2*Pi,  frames= 54,  orientation= [40,70], shading= z, axes= box, style= patchnogrid ;

@Carl Love  I think that  discont  option  can only build one circle at the respective values  x. There is no need to build a hollow circle on the right end, because we consider the single-valued functions. For example, when  x = 2, the function value is 2, and the solid circle indicates it .

@Markiyan Hirnyk  Thanks. I think everyone can come up with many applications of the procedures  Composition  and  Composition1. Here is one of them:

The problem - how many ways 1 dollar (100 cents) can be exchanged  by 10 coins (in denominations of 1, 5, 10, 25, 50 cents) in any combination (exactly 10 coins should be used)?

Composition1(100,10, {1,5,10,25,50}):

convert(map(t->sort(t), %), set);

nops(%);

@Markiyan Hirnyk  for your brilliant solution!

@Markiyan Hirnyk  If my explanation is not enough for you, here's the direct check. For it, I added to the original equation 3 additional equations  diff(y(x), x$2)=k(x), diff(k(x), x)=p(x), diff(p(x), x)=q(x)

restart;

Digits:=20:

z:=x->diff(y(x),x,x)+sin(y(x)):

u:=diff(z(x),x):

v:=diff(z(x),x,x):

w:=diff(y(x),x$3)+diff(y(x),x$2):

sys:=subs({y(x)=a,diff(y(x),x)=b,diff(y(x),x$2)=c,diff(y(x),x$3)=d,diff(y(x),x$4)=0.1},{z(x)=0,u=0,v=0, w=-0.1});

fsolve(sys);

assign(%);

sol:=dsolve({z(x)=0,diff(y(x),x$2)=k(x),diff(k(x),x)=p(x),diff(p(x),x)=q(x),y(0)=a,D(y)(0)=b,k(0)=c,p(0)=d}, numeric);

sol(0);

%[2]+%[3], %[4];

@Markiyan Hirnyk  Let  y(x)  is a solution of the original equation, satisfying to conditions  diff(y(x), x$3)(0)+diff(y(x), x$2)(0)=-0.1  and  diff(y(x), x$4)(0)=0.1

Substitute  y(x) into the original equation and differentiate it 2 times. We got 3 identities that hold for any x, in particular, for x=0. Introduce the notations for  y(0)  ( a=y(0) ) and for all derivatives at  x=0  and adding your initial condition, we obtain a system of 4 equations with 4 unknowns. The further solution is obvious.

You wrote "Is it possible to verify those conditions by the fdiff command?" Directly is impossible, because  y(x)  should be an algebraic expression. A possible approach - using individual points  [x, y(x)]  firstly to get this expression, for example, by interpolation.

@Markiyan Hirnyk  These equations are your initial conditions (to calculate the arbitrary constants   _C1  and  _C2  in symbolic solution  sol  in my first post in this thread).  Also pay attention to the addition to my second post.