I assume that you are loading the PDEtools package because you intend to solve a PDE.  Let mypde be that PDE.  Then instead of solving mypde,  solve the system consisting of { mypde, diff(u,x)=0 }.  This will force diff(u,x) and all its higher order derivatives to be zero throughout.

The short answer to your question is:

q1 := eval(q, omega=2646.408147):
NullSpace(q1);

That said, I should add that I wouldn't put trust in this result.  I fact, you shouldn't put trust in your own matrix q, which contains numerical entries like 10^(-118)  and  10^17.  To put things in perspective, let's note that the distance to the Moon is "only" 3.844 x 10^5 km.  Note the 10^5.  Distance to the Sun is "only" 1.52 x 10^8 km.  Distance to the nearest star is 4 x 10^13 km.  What does your 10^17 measure?  And in what units?

I suggest that instead of continuing with the current calculation (which I consider insane) go back to the derivation of the equations and see where those numbers come from.  If you were doing astronomy, for example, you wouldn't express distance to Neptune in millimeters.  Instead, you would say that the distance to Nepute is 30 AU.  One AU (Astronomical Unit) is the distance between the Earth and the Sun.

Huge variations in numerical coefficients result in numerical instability and nonsensical results.  Choose your units correctly to obtain reliable answers.

Let's see what ZZ looks like as a function of omega:

plot(ZZ, omega=0..10);

We see that the graph confirms your calculation that the positive roots are at:

     1.360526281, 4.077926876, 6.798915756, 9.453815360, 9.864454570

Now, you plug the root 1.360526281 in ZZ and you get something like -1.11*10^58 rather than zero.  There is nothing wrong in that.  Considering that the function's amplitude is of the order of 10^67, a number like 10^58 is one billionth of the amplitude, and therefore it is zero for all practical purposes.

The attached worksheet contains a solution.  I am sure there is room for improvement.

happy-numbers.mw

restart;
with(LinearAlgebra):
A :=<4,1,0,0; 0,0,4,1>;

N := NullSpace(A);

The vectors in N are what you are looking for since they are orthogonal to the A's rows.  Thus we form the larger matrix by appending them to A:
< A, N[1]^+, N[2]^+>;

It's not difficult to automate the operation.  Ask if you don't know how.

To understand the source of the problem, consider the simplified recursion u[i] := u[i-1] + u[i+1].  You begin with i=1 and get:

u[1] := u[0] + u[2]:

No problem thus far, but at the next step where i=2, we get:

u[2] := u[1] + u[3];

In view of the previous definition of u[1],  this expands to

u[2] :=  (u[0] + u[2]) + u[3];

which attempts to define u[2] in terms of itself, leading to an infinite regress.  That's why Maple complains about recursive assignment.

I don't understand why you disallow parametrization, but if you insist, then here is a solution without parametrization:

restart;

with(plots): with(plottools):

# produce a standard circle
c0 := circle([0,0], 1):

# embed it in the x-z plane:
c1 := transform((x,y)->[x,0,y])(c0):

# rotate it about the z axis by Pi/3
c2 := rotate(c1, Pi/3, [[0,0,0], [0,0,1]]):

display([
  sphere([0,0,0], 1),
  display(c1,color=red, thickness=3),
  display(c2,color=blue, thickness=3)
], scaling=constrained, labels=[x,y,z], orientation=[-130,72]);


If, however, you don't object to parametrization, there is a simpler solution:

restart;

with(plots): with(plottools):

# orthonormal vectors in the circle's plane:
u := <cos(Pi/3),sin(Pi/3),0>:  v := <0,0,1>:

# plot the circle
c := spacecurve(u*cos(t) + v*sin(t), t=0..2*Pi):

display([
    sphere([0,0,0], 1),
    display(c, color=blue, thickness=3)
], scaling=constrained, labels=[x,y,z], orientation=[-130,72]);

The trick in producing the oscillating pendulum animation is to set things up so that the system's overall motion is periodic.  If the periods of the individual pendulums are integers, then the period of the overall system equals the least common multiple of those periods.  If the periods are selected so that their least common multiple is not too large compared to each one of them, then the result is a pleasing animation.  The attached worksheet pendulums.mw shows how.  Here is a sample animation.  The motion was calculated over one period,  Since the overall motion is periodic by design, then playing the motion in a loop produces the correct result for all times.

Edit: Updated the earlier worksheet and animation by changing the line

periods := numtheory[divisors](8!)[83..94];

to

periods := numtheory[divisors](13!)[1571..1582];

to produce an even nicer animation.

I looked through your worksheet but did not attempt to debug (I find Maple's 2D input hard to read.)  Instead, I just implemented the book's algorithm in Maple myself. Its result agrees with the book's.

burden-faires-algo-12.1.mw

As I commented earlier, inertial forces are essential in forming Lagrangian points.  Let's say we have a two-body system of masses a and 1-a, rotating in circles about their common center of mass. In a Cartesian coordinate system that rotates with the system, the masses are located at (a, 0) and (1 - a, 0), therefore their combined gravitational potential is

( 1 - a ) / sqrt( ( x - a )^2 + y^2) + a / sqrt( ( x + 1 - a )^2 + y^2 ) ;

A particle of neglibible mass placed within this rotating system will be subject to the gravitational pulls of these two masses as well as a centripetal force that grows linearly with distance from the origin.  The potential of that force will be a quadratic, therefore the particle will be subject to the effective potential field

V := ( 1 - a ) / sqrt( ( x - a )^2 + y^2) + a / sqrt( ( x + 1 - a )^2 + y^2 )  + ( x^2 + y^2 ) / 2 ;

Now let's set a=0.1 and plot the level curves of V:

plots:-contourplot(eval(V,a=0.1), x=-2..2, y=-2..2,
                               contours=[seq(i, i=1..3, 0.1)], grid=[100,100]);

In the definition of myeq, change r and s to r(y) and s(x).

Aside: I haven't tested this because what you have posted are picture images which I can't copy and paste into a Maple worksheet.  It would have been more helpful if you would have posted the textual forms of your inputs for those expressions.

That integral as stated is undefined:

• Note the -lambda*t+1 in the denominator.  What happens when t hits 1/lambda?
• What happens to exp(-beta*t) as t goes to infinity?  It depends on the sign of beta!
• What happens to t^a as t goes to zero?  It depends on the sign of a!

You need to make assumption to get out of that mess.  Your integrand is:

So the integral portion of your expression evaluates to:

int(t^a*exp(-beta*t)/(-lambda*t+1), t = 0 .. infinity) assuming beta>0, lambda<0, a > -1;

Aside: The Gamma function is spelled GAMMA (all capitals) in Maple.

collect(B, [x,y], distributed, factor);

or perhaps

collect(B/2, [x,y], distributed, factor);

I don't know a rainbow's true color distribution.  Here is a fake:

plot3d(3-r, t=0..Pi, r=1..2, coords=cylindrical,
    scaling=constrained, style=surface, shading=zhue,
    lightmodel=none, axes=none, orientation=[-90,0,0]);