In the attached worksheet I show how to generalize the definition of the pedal curve to pedal surface and calculate the equation of the pedal surface in the general case.  Then I illustrate the result by calculating the pedal surface of an ellipsoid with respect to its center.  Here is what it looks like.

Download pedal-surface.mw

No need for extra packages.  Just do this:

restart;
k := x -> piecewise(x < -6, 1/2*x + 10, x < -3, 1, x < 6*Pi, sin(1/2*x) + 2, x < 25, x - 15);
plot3d(k(x), t=0..2*Pi, x=-20..20, coords=cylindrical);

The parametrization of a surface (any surface) is not unque.  The way a surface is parametrized can have a drastic effect in how it is rendered on screen.

The standard parametrization of the oloid does not produce good computer-generated graphics as you have observed. But an oloid has a great deal of symmetries.  It is not difficult to see that it may be decomposed into four congruent patches, each of which has a quite "tame" parametrization for which Maple's default plotting grid is more than adequate.  Here I construct the oloid by putting together the four patches. 

restart;
alpha[1] := sin(t):
alpha[2] := -1/2 - cos(t):
alpha[3] := 0:
beta[1] := 0:
beta[2] := 1/2 - cos(t)/(1+cos(t)):
beta[3] := sqrt(1 + 2*cos(t))/(1+cos(t)):
x := (1-m)*alpha[1] + m*beta[1];
y := (1-m)*alpha[2] + m*beta[2];
z := (1-m)*alpha[3] + m*beta[3];
plot3d([[x,y,z], [x,y,-z], [z,-y,x], [-z,-y,x]], m=0..1, t=-Pi/2..Pi/2,
    scaling=constrained, color=[red,green,blue,yellow]);

I have colored the patches in four different colors in order to delineate them.  Change surface style and color as it suits your needs.

T := t -> (t+459.47) * 5/9;
seq(printf("%10d %10.2f\n", t, T(t)), t=100..110);
       100     310.82
       101     311.37
       102     311.93
       103     312.48
       104     313.04
       105     313.59
       106     314.15
       107     314.71
       108     315.26
       109     315.82
       110     316.37

Inversion in a sphere is much simpler than what you have in your worksheet.  Here is how I would do it.

Download inversion-in-sphere.mw

There are mnay ways of doing it, depending on what you want to do with the results.  Also it depends on how you represent quaternions.  You say that you want to represent a quaternion as a vector of length 4, but in your example your write qaa and qab as lists.  Lists and vectors are different types of objects in Maple.

Let's say you represent a quaternion as a list of length 4.  Then this will do what you want:

That said, I much prefer the [scalar, vector] pair representation of quaternions (see the Wikipedia article) rather than lists of length 4.  But that again depends on what you want to do with them.

de := a*diff(y(x),x,x) + b*diff(y(x),x) + c*y(x) = 0;
dsolve(de) assuming b^2 < 4*a*c;

Download mw.mw

Okay, your hand-written calculations are good.  Here I will do those in Maple.

Download mw.mw

How do you like this?

restart;

# macro emulates Matlab's "jet" colorscheme:
macro(jet=["zgradient", ["DarkBlue", "Blue", "DeepSkyBlue", "Cyan", "LawnGreen", "Yellow", "Orange", "Red", "DarkRed"]]):

plot3d((1-xi^2)*(1+tau), tau=-1..1, xi=-1..1,
    colorscheme=jet, orientation=[-120,72,0], axes=framed);

plot3d((1-xi^2)*(1+tau), tau=-1..1, xi=-1..1,
    colorscheme=jet, gridstyle=triangular, orientation=[-120,72,0]);

As Kitonum has noted, this is not an infinite fraction.  The integrand is z = ..  Then we see that z^2 = x*z, and therefore z = x.  Then the integral of z is x^2/2 + constant.  QED.

Forget about solve.  After you define psi, just do:

plots:-contourplot(psi, x=-1..1, eta=-2..2);

Try this:

restart;
plot(x^2, x=-1..1, labels=[typeset(NFE), typeset(log__10(`Max Err`))]);

Here is a slightly modified version of your worksheet.  The determinant grows too fast for efficient numerical calculations.  Therefore I have scaled it by dividing by cosh(...) which does not affect the roots since cosh is never zero.  Then you see that you can calculate the first five roots without any fuss.

Download Natural_frequency_No_Foundation_Mass-modified.mw

If you look closely at your equations, you will see that at some places you have x, and at some other places you have x(t).  They should all be x(t).  Same goes for y(t) and u(t).

I have fixed these in the attached worksheet.  Note the technique—to define a differential equation dx/dt = F(x,y,u), we begin with defining the function (not an expression!) F, and then express the equation as dx(t)/dt = F(x(t), y(t), u(t)).

But note that the solution to your system is not defined in the entire length of the time interval that you have requested.  Since you have x, y, and u in the denominators, the solution stops when any of these variables reaches zero.

Download mw.mw