Here I do Exercise #7.  You do the rest.

Download mw.mw

This is a linear system, therefore (0,0) is an equilibrium point.

To plot a phase portrait, do:

DEplot(sys, [x(t),y(t)], t=0..1, x=-4..4, y=-4..4, arrows=comet, color="DarkCyan",
    dirfield=[15,15], size=[800,600], axes=boxed);

I much prefer quaternions over Euler angles for expressing rotation of rigid bodies.  Euler angles work well when the rigid body's orientation is limited to rather narrow range (think of a spinning top).  Quaternions work well when a rigid body has the possibility of turning in every which way (think of a satellite in space).

Stylistically, each of the three Euler angles  serves a different purpose.  In contrast, the components of a quaternion participate symmetrically in the equations of dynamics, resulting in a more pleasant treatment.

That said, Euler angles tend to result in a more compact form of the differential equations of motion.  These are preferable for "calculations by hand".  The differential equations of motion obtained via quaternions are generally more complex, and are better suited for numerical calculations on a computer.

Here is a mini-lession on quaternion algebra

A quaternion is a pair [a, u], where a is a scalar and u is a 3D vector.  The algebra of quaternions prescribes the algebraic properties of such a pair.  Thus, the sum of the quaternions [a,u] and [b,v] is defined in the expected way as [a+b, u+v].  The scalar product  of a number c and a quaternion [a,u] is [ca, cu].

The interesting new concepts are those of the conjugate [a,u]* of a quaternion and the product [a,u] o [b,v] of two quaternions, defined as

[a, u]* = [a, -u]
[a, u] o [b,v] = [ab - u.v, av + bu + uxv]

where u.v and uxv are the dot product and the cross product of the vectors u and v.  Note that the quaternion product is not commutative because uxv = - vxu.

That's pretty much all of the quaternion algebra needed for expressing rotations in 3D.

To see that, let n be an arbitrary unit vector in 3D and phi an arbitrary angle.  Define the rotation quaternion

q = [ cos(phi/2), n sin(phi/2)]

Let P be an arbitrary point P in space and let r be its position vector, that is, the vector that stretches from the origin to P.  Then form the quaternion product

[a,v] = q o [0,r] o q*

Expanding the product according to the quaternion rules laid out above, we find that a = 0, that is

[0,v] = q o [0,r] o q*

Furthermore, an elementary geometric analysis shows that  the vector v is the result of rotating the vector r about the vector n by the angle phi.

Thus, to rotate a rigid body about a vector n by angle phi, apply the formula obtained above to all of its points. That expresses the rigid body's rotation in a natural way; there are no Euler angles!

End of the lesson

To connect this to what you have shown in your post, MapleSim expresses the rotation quaternion q as a list of four numbers, q0 for the scalar part, and q1, q2, q3 for the components of the vector part.  Since n is a unit vector in the rotation quaternion, we have q0^2 + q1^2 + q2^2 + q3^2 = 1.  That's the "algebraic constraint" seen in MapleSim's message.  There is some confusion there about the rest of the message.  There are 7 (not 8) generalized coordinates, x, y, z, q0, q1, q2, q3.  The duplicate q0 at the end of the vector vQ should not be there, unless MapleSim has some good reason for expressing vQ that way.

I can only guess what you are asking.  Consider formulating the question in your native language and having Google Translate change it to English.

Here is my best attempt to make something meaningful from what you have written.

The function f given in your post has infinitely many roots.  The roots are grouped in clusters.  You wish to find the first cluster which contains 50 or more roots.

It turns out that the first cluster contains 122 roots, so it meets your requirement.

To see that, we change the variable from x to s by setting x^3 = s that is, x = s^(1/3).  Then the roots (in s) within each cluster will be more or less uniformly distributed, making it easier for fsolve() to find them.  Here are the details.

Download mw.mw

I can't tell why Maple does what it does.  But a close look at the Maple's result makes it possible to find the correct solution as follows.

Download mw.mw

Maple's mod, modp, and mods are intended for use in polynomial algebras.  For modular arithmetic with integers use irem instead.  Thus,

FB := sum(irem(5, product(2, t = 0 .. irem(q - 1, 3))), q = 1 .. 2);

evaluates to 2 as you would expect.

>  restart; >  with(plots): >  HarmonicEnergies := proc(N::posint)     local n;     return Vector([seq(n + 1/2, n = 0..N-1)]); end proc: >  BuildM := proc(N::posint, mx::positive)     local m, n, Mentry, M;     Mentry := (m,n) -> exp(1/(4*mx)) *         sqrt(min(m,n)!/max(m,n)!) *         (1/sqrt(2*mx))^abs(m-n) *         LaguerreL(min(m,n), abs(m-n), -1/(2*mx));     M := Matrix(N, N, (m,n) -> evalf(Mentry(m-1, n-1)));     return M; end proc: >  GaussianStateCoeffs := proc(N::posint, mx::positive,                              x0::realcons := 0.0, vx0::realcons := 0.0,                              omega::realcons := 1.0,                              sigmax::realcons := 1/sqrt(2*mx*omega))     local px0, alpha, pref, C, n;     px0 := mx * vx0;     alpha := x0 / (sqrt(2) * sigmax) + I * px0 * sigmax;     pref := exp(-abs(alpha)^2 / 2);     # Proper n from 0 to N-1     C := Vector(N, n -> evalf(pref * alpha^(n-1) / sqrt((n-1)!)), datatype = complex[8]):     return C; end proc: >  N := 3: mx := 2: my := 1: lambda := 1: y0 := 0: p0 := -5: t_final := 1:   # was 30(!) >  E := HarmonicEnergies(N); >  M := BuildM(N, mx); >  B := Matrix(N,N, (m,n) -> exp(I*(E[m]-E[n])*tau) * M[m,n]); >  Cvec := Vector(N, i -> C[i](tau)); >  de_C := seq(diff(Cvec,tau) =~ lambda*exp(-y(tau)) * B . Cvec); >  de_y := diff(y(tau),tau) = p(tau) / my; >  de_p := diff(p(tau),tau) = lambda * exp(-y(tau)) * Cvec^+ . B . Cvec; >  ic := seq(subs(tau=0, Cvec) =~ GaussianStateCoeffs(N, mx)),         y(0)=y0, p(0)=p0; >  dsol := dsolve({de_C, de_y, de_p, ic}, numeric); Plot the real and imaginary parts of y(tau) >  display(     < odeplot(dsol, [tau, Re(y(tau))], tau=0..t_final) |       odeplot(dsol, [tau, Im(y(tau))], tau=0..t_final) > );         >    >  display(     < odeplot(dsol, [tau, Re(p(tau))], tau=0..t_final) |       odeplot(dsol, [tau, Im(p(tau))], tau=0..t_final) > );         >    >  display(     < odeplot(dsol, [tau, Re(C[1](tau))], tau=0..t_final) |       odeplot(dsol, [tau, Im(C[1](tau))], tau=0..t_final) > );         >    >  display(     < odeplot(dsol, [tau, Re(C[2](tau))], tau=0..t_final) |       odeplot(dsol, [tau, Im(C[2](tau))], tau=0..t_final) > );         >    >  display(     < odeplot(dsol, [tau, Re(C[3](tau))], tau=0..t_final) |       odeplot(dsol, [tau, Im(C[3](tau))], tau=0..t_final) > );         >

Download mw2.mw

I assume that you are writing your notes in LaTeX.  If not, then you should.

You may do the drawings in Maple, but that will be a struggle.  Things will go much smoother if you learn how to use the TikZ package and do the drawings in LaTeX.  That will require some investment of time and effort, but it would be a worthwhile investment.

Here is a sample.  Took me less than 10 minutes to do it.

\documentclass{article}
\usepackage{tikz}
\begin{document}

\begin{tikzpicture}[>=latex]
    \pgfmathsetmacro{\R}{2.5}
    \pgfmathsetmacro{\h}{0.05*\R}

    % the coordinate axes
    \draw[->] (-1.2*\R,0) -- (1.2*\R,0) node[below] {$x$};
    \draw[->] (0,-0.1*\R) -- (0,1.2*\R) node[right] {$y$};

    % the contour
    \draw[red, thick]
        (-\R,0) 
        -- (-1/\R,0)
        arc [start angle = 180, delta angle = -180, radius = 1/\R]
        -- (\R,0)
        arc [start angle = 0, delta angle = 180, radius = \R]
        -- cycle;

    % arrowheads
    \draw[->, >=stealth, red, very thick] (-0.5*\R,0) -- +(0.01,0);
    \draw[->, >=stealth, red, very thick] ( 0.5*\R,0) -- +(0.01,0);
    \draw[->, >=stealth, red, very thick] ( 0, \R)    -- +(-0.01,0);

    % labels
    \node at (45:1.15*\R) {$\mu_R$};
    \draw
        (-\R,0) -- +(0,-\h) node [below] {$-R$}
        (-1/\R,0) -- +(0,-\h) node [below] {$-\frac{1}{R}$}
        ( 1/\R,0) -- +(0,-\h) node [below] {$ \frac{1}{R}$}
        ( \R,0) -- +(0,-\h) node [below] {$ R$}
    ;
        
\end{tikzpicture}
\end{document}

>  restart; >  with(plots): The vertices: >  A, B, C := <0,0>, <3,0>, <3*cos(Pi/3), 3*sin(Pi/3) >; The points P, Q, R: >  Q := 1/3*A + 2/3*C; P := 1/3*C + 2/3*B; R := <4,0>; If P lies on QR, them this equation will have a solution, otherwise there won't be a solution: >  P = (1-t)*Q + t*R; solve(%, t); >  display(         pointplot([A,B,C,A], connect, color="Green"),         pointplot([B,R], connect, color="Green", linestyle=dash),         pointplot([P,Q,R], symbol=solidcircle, symbolsize=20, color="Red"),         pointplot([Q, R], connect, color="Blue"),         textplot([                 [seq(A), 'A', align={below,left}],                 [seq(B), 'B', align={below}],                 [seq(C), 'C', align={above}],                 [seq(Q), 'Q', align={above,left}],                 [seq(P), 'P', align={above,right}],                 [seq(R), 'R', align={below}],         NULL]), scaling=constrained, axes=none);

Download mw.mw

Download mw.mw

Added later:

With the proc(m,n,col) on hand, we can generate tilings of desired size:

M, N := 10, 7:
seq(seq(p(m,n,red),   n=irem(m,2)..N, 3), m=0..M),
seq(seq(p(m,n,green), n=irem(m,2)+1..N, 3), m=0..M),
seq(seq(p(m,n,blue),  n=2*irem(m+1,2)..N, 3), m=0..M):
display(%, scaling=constrained, axes=none);

Here is yet another solution.  The logic is explained in the attached PDF.

area-puzzle.pdf

Download area-puzzle.mw

The "solution" that you have obtained says that y(t)^2 = arctan(t) / t.   Letting t go to zero we get y(0)^2=1.  In what way does that satisfy the initial condition y(0)=0?

If you are looking for a way to mark a missing point on a graph by an empty circle, as it's done in textbooks, this may be of some use to you.

>  restart; >  with(plots): >  expr := (x^3-1)/(x-1);    # change as needed Want to remove the point (x0,y0) from the graph >  x0 := 1; >  y0 := limit(expr, x=x0); Size of the gap in the graph >  h := 0.025; The symbolsize is determined by trial-and-error >  display(   plot(expr, x=0..x0-h),   plot(expr, x=x0+h..2),   pointplot([x0,y0], symbol=circle, symbolsize=30), color="Green", thickness=3);

Download mw.mw

You have been shown several fancy approaches to solving your problem. Here is a  very basic approach, which is not much different from what you already have in your worksheet.  It's just organized in a cleaner way.

Download mw.mw

Consider using the inert form %* of multiplication: