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

I think you are asking for a surface of revolution.  A tubeplot is something entirely different.

A surface of revolution may be ploted in Maple in a very many different ways, the simplest of which is with the help of Student[Calculus1] package, like this:

restart;
with(Student[Calculus1]):
SurfaceOfRevolution(x^2+x-1, x=0..1, axis=vertical, output=plot, axes=boxed);  p1 := %:

We see that the vertical axis goes up to z=1.  You want to continue this straight up, to height z=5.  The graph of that extra part is an ordinary cylinder, which we produce through

plot3d(1, t=0..2*Pi, z=1..5, coords=cylindrical, color="Maroon", style=surface);  p2 := %:

Now p1 and p2 represent thos two plots. We put them together with the help of plots:-display command:

plots:-display([p1,p2], scaling=constrained);

For some reason which I don't understand, this does not do the expected thing.  It shows only the p1 graph  This may be a bug in Maple, or perhaps something else that I am doing wrong.  I will let others comment on that.

Added later:

I just noticed that if we specify a view=... option to the latter command, it displays the composite picture properly:

plots:-display([p1,p2], scaling=constrained, view=-1..6);

The reason that it does not work without the view is that the SurfaceOfRevolution command sets the view option to its own extents, which is view=-1..1 in this case.  That option caries over to the plots:-display command, and therefore the result is that the composite surface is truncated at the z=1 level.  We specify the view=-1..6 option explicitly to override the one set by SurfaceOfRevolution.

The way I see it, this is sort of a minor bug, but perhaps not worth making a big fuss about.

You don't need DEtools for this, but you need plots:

restart;
with(plots):
w := 1+t/100;
sys := {diff(x(t), t) = y(t), diff(y(t), t) = -w^2*x(t), x(0) = -1, y(0) = 2};
dsn := dsolve(sys, numeric);
odeplot(dsn, [[t,x(t)],[t,y(t)]], t=0..10);
E := 1/2*(y(t)^2 + w^2*x(t)^2);
odeplot(dsn, [t, E], t=0..10);
J := E/w;  # I is predefined in Maple; use J
odeplot(dsn, [t, J], t=0..10);

I have never felt the need to use animate() whose syntax I find rather cluncky.  I just prepare the sequence of the individual frames first, then use display() to animate it, as in:

with(plots):

frames := seq(
  dualaxisplot(plot(sin(x+A), x=0..5), plot(cos(x+A), x=0..5)),
  A=0..5
):

display([frames], insequence=true);

I couldn't execute your worksheet in Maple 2015.  It ran into errors.

But just looking through it I formed an idea of what it is that you want to do and produced this worksheet: mw.mw

The end result is:

As to your question about the area of the spherical triangle, there are very many neat formulas for it.  See the Wikipedia article:

    https://en.wikipedia.org/wiki/Spherical_trigonometry

where you wil find formulas for the area in terms of the angles a,b,c (these are the angle between the vectors shown in the graphics above) or the angles A,B,C (these are the actual angles at the triangle's vertices measured on the sphere's surface).  Both sets of angles may be calculated easily in terms of dot products and cross products from the given data.

If I am not mistaken, the following worksheet does what you are asking.

mw.mw

Here is the graph of the function which is expressed as an integral in your post:

Here is one of many different ways of writing a loop:

for V from 1 to 30 do
   ....
end do;

convert(abs(1-b)+abs(1+b), piecewise);

D[1,1](w) means the second derivative of w(x,y) with respect to x.  Then D[1,1](w)(x,0) means the second derivative in the x direction along the bottom edge.  Surely you don't mean that.  If your intent is to specify zero bending moment along the bottom edge, you want D[2,2](w)(x,0)=0.

General advice: If you are having difficulty with a complex problem, see if you can do a simpler problem first.

So forget about your F_T for the moment.  Do you know how to work with your Esolve()?

For that, examine what Esolve() does, and what its parameters are.

We see that Esolve() is designed to solve a system of differential equations of the form dy/dx = f(y,x) for the unknown vector-valued function y(x).

The intial condition is y(x0) = y0.

The solution is produced by taking N steps of stepsize h.

You are interested in solving a problem of the type dv/dt = F(v,t).  Thus, your v and t correspond to Esolve()'s y and x, respectively.

Let's give it a try with a function f which takes a vector <v[1], v[2]> as argument and returns the vector <-v[2],v[1]>.  Thus, in Maple we set:

f := v -> <-v[2], v[1]>;

For the initial value pick v(0) = <1,0>.  Now solve with stepsize h = 0.1 and take 25 steps:

sol := Esolve(f, 0.1, 0, <1,0>, 25);

Plot the solution, that is, v(t) versus t:

plots:-pointplot([seq(sol[i], i=0..25)], scaling=constrained, color=red, connect=true);

See if you can take it up from there.

This solution is obtained by eyeballing the figure, not through any algorithm: