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

The wording of your question is a bit ambigious, so let's see if I can clarify.

Let S be the surface in R^N defined by A*x = b, where A and b are as you have stated.  Let p be a point in R^N.  We wish to determine the point q which is the orthogonal projection of p onto the surface S.

Without loss of generality you may assume that rank(A) = d, otherwise do a row echelon reduction of the system Ax=b, and then discard the entirely zero rows, if any.

It can be shown that:

q = [ I − A^T (A A^T)⁻¹ A ] p + A^T (A A^T)⁻¹ b.

where I is the N × N identity matrix, and the T superscript indicates the transpose.  I needed such a formula a few years ago as a part of the book that I was writing.  I couldn't find a proof in a brief literature search, so I derived it myself.  I wouldn't be surprised if it turns out that this is an old and well-known formula.  In fact, I would appreciate it if someone points me to a reference.

In the mean time if you need a proof, it is on page 213 of my book Programming Projects in C.

See this worksheet:

mw.mw

Change the line

odeplot( sol,[x, y(x)], x=a..b);  # I plot the solution in the interval (a,b)

to

P0 := odeplot( sol,[x, y(x)], x=a..b);  # I plot the solution in the interval (a,b)

Then do

plots:-display(P0,P1,P2);

You have terminated almost every command in your worksheet with a colon, in effect depriving yourself from the benefit of Maple's feedback.

Change the colons to semicolons to see what Maple thinks of your input.  When you do that, you will see that the fourth equatiuon in Eqns is not what you would have expected.  Examining the cause, you will find out that you have entered vyH where you should have vyH(t).

Carl has pointed out several syntax errors and other issues with your equations.

I believe there are more problems.  If the equations are to model an oscillator, the x on the right-hand side of the first equation should be z, and the y on the right-hand side of the second equations should be x.  So we change the right-hand sides to:

f := (x,z) -> alpha*(1 - (x^2/a^2 + z^2/b^2))*z + w*a/b*z;
g := (x,z) -> beta*(1-(x^2/a^2 + z^2/b^2))*x + w*b/a*x;

Then the differential equations are:

de1 := diff(x(t),t) = f(x(t),z(t));
de2 := diff(z(t),t) = g(x(t), z(t));

We introduce the parametes

params := alpha=1, beta=1, a=0.4, b=0.2, w=1;

and evaluate the equations accordingly:

sys := eval([de1, de2], [params]);

To get an idea of the region of interest in the x-z plane, we find the system's equilibria:

RealDomain[solve](
    eval([f(x,z)=0, g(x,z)=0], [params])
);

We see that in addition to the origin (0,0), there are four equilibria at x=0, z=+/- 0.35, and x=+/-0.49, z=0.  Thus we set

xmax := 0.6;  zmax := 0.5; tmax := 8;

and generate a reasonable set of initial data:

 ic1 := seq([x(0)=h,z(0)=0], h=-xmax..xmax, 0.1),
      seq([x(0)=h, z(0)=-zmax], h=-xmax..0, 0.05),
      seq([x(0)=h, z(0)=zmax], h=0..xmax, 0.05);

Now we call DEplot to produce a phase portrait:

DEtools[DEplot](sys, [x(t),z(t)], t= 0..tmax, [ic1],
  linecolor=black, thickness=1, stepsize=tmax/(3*48),
  x(t)=-xmax..xmax, z(t)=-zmax..zmax, arrows=none);

Lo and behold, the Starbucks siren emerges!

We may want to add extra orbits to cover the chest area:

ic2 := seq([x(0)=h,z(0)=0], h=-xmax..xmax, 0.1),
      seq([x(0)=h, z(0)=-zmax], h=-xmax..0, 0.05),
      seq([x(0)=h, z(0)=zmax], h=0..xmax, 0.05),
      seq([x(0)=0, z(0)=h], h=-0.34..0.34, 0.04);

DEtools[DEplot](sys, [x(t),z(t)], t= 0..tmax, [ic2],
  linecolor=black, thickness=1, stepsize=tmax/(3*48),
  x(t)=-xmax..xmax, z(t)=-zmax..zmax, arrows=none);