@dharr has shown how to implement your textbook's approach to solving the PDE.  I just want to remark that a solution may be obtained directy in Maple through the traditional apprach.  Maple's solution appears to be computationally more efficient since it involves an integration over the interval [0,t] as opposed to the textbook's solution which integrates over [0,infinity].

Plotting the result at the default 25x25 grid size takes about 10 minutes on my computer.  On a 10x10 grid it takes about 20 seconds.

Download heat-eq2.mw

Here is a general proc for plotting arrowheads. 

>  restart; >  with(plots): Proc arrowhead()         Draws an arrowhead at the point p in Cartesian coordinates.         The arrowhead's reference point is its reentrant vertex. Arguments:         p: is the arrowhead's location specified as a list   or a vector  .         angle: is the angle of the arrowhead (rotation) relative to the horizontal axis         mag: (optional, default=1) is the magnification factor.  With the default ,         the distant between the reentrant corner and the arrow tip is 1/10.         Any additional arguments are assumed to be plotting options and are passed         to plots:-display. >  arrowhead := proc(p::{[realcons,realcons],'Vector'(2,realcons)}, angle, {mag:=1})         uses plots, plottools;         polygon([[0,0], [-1/2,-sqrt(3)/2], [2,0], [-1/2,sqrt(3)/2]],                 style=polygon);         rotate(%, angle);         scale(%, mag/20, mag/20);         translate(%, p[1], p[2]);         display(%, _rest); end proc:       Example >  arrowhead(<3,4>, Pi/3, mag=10, color=red, scaling=constrained); Example >  A, B, C := <-1,0>, <0,0>, <1,0>; (1) >  display(         pointplot([A, B, C], symbol=solidcircle, symbolsize=25),         plottools:-line([-2,0],[2,0]),         arrowhead((A+B)/2,   Pi, mag=1.2),         arrowhead((B+C)/2,   0,  mag=1.2),         arrowhead(A-<1/2,0>, 0,  mag=1.2),         arrowhead(C+<1/2,0>, Pi, mag=1.2), scaling=constrained, tickmarks=[0,0], size=[800,100]);

 

Download arrowhead.mw

The terms "chaos" and "chaotic" have precise technical meanings.  The equations in your textbook may have nonperiodic solutions, but they are not chaotic.  Furthermore, the textbook refers to the u--v diagrams as "phase portraits".  But they are not phase portraits.  There is no meaning to a "phase portrait" for nonautonomous systems.  [Nonautonomous means that the independent variable (eta in your case) appears explicitly on the right-hand sides of the equations.]  With that in mind, here are a few things that you can do with that system.

Download mw.mw

Why Mathematica or Matlab?  This is a Maple forum, so you will get a Maple answer.

The easiest way to plot the system's phase portrait—essentially just two lines of Maple code—is to observe that your system of differential equations is conservative and its potential is

P := 1/2*F1*u^2 + 1/3*F2*u^3 - 1/2*v^2;

The system's orbits are the level curves of P, so here is what we get:

subs(F1=1, F2=-1, P):
plots:-contourplot(%, u=-1.5..2, v=-1.5..1.5, scaling=constrained,
    colorscheme="DivergeRainbow", contours=[seq](x, x=-0.4..0.4,0.1));

In a reply to your previous question, I explained how to upload your worksheet.  Don't be shy.  Uploading your worksheet will help others to help you.  Anyway, here is a solution, shot in the dark.  

Give a name to your animation, as in 

QQ := animate(...);

Then do:

with(DocumentTools): with(DocumentTools:-Layout):
InsertContent(Worksheet(Group(Textfield(InlinePlot(QQ, scale=2.5)))));

Change the 2.5 to an appropirate scale factor.

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