@acer   Oh sure. My goal was to keep the idea of initial code, but to make it work.

@Zeineb I think it's possible, but I do not know how, because I do not work with Matlab.

A quote from the help:

"Calling Sequence
ProbabilityFunction(X, t, options)

Parameters
X - algebraic; random variable or distribution
t - algebraic; point (assumed to be an integer) ...."

@Zeineb 

Theta:=[ 0, Pi/3, 2*Pi/3, Pi]:
Phi:=[ 0, Pi/3, 2*Pi/3, Pi, 4*Pi/3, 5*Pi/3, 2*Pi, 7*Pi/3]:
F:=(theta,phi)->[cos(phi)*sin(theta),sin(phi)*sin(theta),cos(theta)]:
{seq(seq(F(theta,phi), phi=Phi), theta=Theta)}:
A:=plots[pointplot3d](%, color=red, symbol=solidcircle, symbolsize=15):
B:=plottools[sphere]([0,0,0], 1, style=surface):
C:=plots[spacecurve]([cos(t),0,sin(t)], t=0..2*Pi, color=blue):
C1:=plots[spacecurve]([cos(t)*cos(Pi/6),sin(t)*cos(Pi/6),1/2], t=0..2*Pi, color=blue):
C2:=plots[spacecurve]([cos(t)*cos(Pi/6),sin(t)*cos(Pi/6),-1/2], t=0..2*Pi, color=blue):
plots[display](A, B, seq(plottools[rotate](C,Pi/3*k,[[0,0,0],[0,0,1]]), k=1..6), C1, C2, axes=normal, view=[-1.4..1.4,-1.4..1.4,-1.4..1.4]);
                

Edit.

@Markiyan Hirnyk 

B2 = (1/2)*sqrt(27*z^3-27*z^2+4) - (3/2)*z*sqrt(3*z-3) = 

= sqrt((27*z^3-27*z^2+4)/4) - sqrt((27*z^3-27*z^2)/4) >0  because  sqrt  is an

increasing function and   27*z^3 - 27*z^2>0  if  z>1

@Markiyan Hirnyk  I don't understand what your (-1)*(-1);  means ?

Under the condition  z>1  it's obvious that  B1>0  and  B2>0

@Markiyan Hirnyk   If  B1*B2 = 1  (see above) then  B1^(1/3) * B2^(1/3) = 1

@Lali_miani  I've already wrote  "To simplify the answer the numerical values of the system parameters should be specified" 

For example:

indets(Sys) minus {x, y, z};

L :=  %=~ [1, 2, 3, 4, 5, 6, 7, 8, 9];

solve(eval(Sys, L), {x, y, z});

evalf([%]);

@tomleslie  When I simplified the specification of the system, I accidentally deleted this factor. I've fixed it.

@one man   The package  geom3d  does not contain any direct commands for plotting and finding the parameters of a circle in 3D, so this work is very useful. The value of the work will be higher if it to be written as a procedure for any three points in 3D (of course, if these points do not lie on a straight line).

@Spinosaurus  I do not see any problems. You yourself put tickmarks on  y-axis  which you need for a corresponding plot:

P:=plot([sin(x)/4, cos(x)+2.4, 1.2], x=0..2*Pi, -1.2..3.6, color=[red,blue,black], tickmarks = [piticks, [seq(0.25*i = i*10^(`-5`), i = -1 .. 1), seq(0.5*i+2.4 = 0.5*i, i = -2 .. 2)]], axes=box):
T:=plots:-textplot([[4.2,0.2,y=sin(x)/10^`5`], [4,2.5,y=cos(x)]], font=[times,roman,14]):
plots:-display(P, T, scaling=constrained);

@asa12   Add the option  linestyle = solid  to the code.

@asa12  I repeat again:  you must first click by the mouse on the plot, and press "play " (a triangle symbol) on the animation panel.

@asa12  In my code  sols  is  a list not a vector. So try first convert your  sols[11]  into list:

convert(sols[11], list);

@asa12   My code works in Maple 12. You should just click by the mouse on the plot, and press "play " on the animation panel.
n is the number of vertices in the broken line (nops(sols)= 9), and 5 is the number of frames to display each segment. If the number of frames to increase, it will be more smooth animation.