This seems to give a very good fit for your data.

> F:= x -> 
      a[0] + a[1]*(x-c)^d + a[2]*(x-c)^(2*d) + a[3]*(x-c)^(3*d);

> Optimization[LSSolve]([seq(F(ET[i,1])-ET[i,2],i=1..21)],
   c=-8 .. -6.25629, d = 0 .. 1, a[0]=-300..-200,
   a[1] = -200..200, a[2] = -2..2, a[3] = -1..1);

[.273160086892695775, [c = -6.25630480630643, d = .606531397333057, a[0] = -269.045867172506, a[1] = 90.3343384288194, a[2] = -1.02769900432745, a[3] = .140313894677253]]

> F0:= subs(%[2],eval(F));

F0 := proc (x) options operator, arrow; -269.04586717250555+90.33433842881936*(x+6.256304806306433)^.6065313973330575+(-1)*1.0276990043274528*(x+6.256304806306433)^1.213062794666115+.14031389467725272*(x+6.256304806306433)^1.8195941919991725 end proc

> with(plots):
  display([pointplot(ET, symbol=cross), 
      plot(F0(x), x = -6.256305 .. 18.3)]);

> residuals:= [seq(ET[i,2] - F0(ET[i,1]), i = 1 .. 21)];

residuals := [0.455784686386096e-2, -.409871157119142, 0.483430078498088e-1, .257717667691310, .340925329193453, 0.587866928924825e-1, -0.964362112382560e-1, -.202394168662721, -.182473165816958, .259676555775968, .115848253096310, -0.877973834528447e-1, -0.786415061487276e-1, -0.283169341274743e-1, -0.397643871756941e-1, -0.290870868733748e-1, 0.358773900078972e-1, 0.484861675207071e-1, 0.275048513215097e-1, 0.173396569728084e-1, -0.602814185723446e-1]

> DynamicSystems[DiscretePlot](ET[...,1], residuals, style = stem);

It's only going through it once because you only told it to go through once.  You need a nested loop, something like:

for j from ... to ... do
  for i from ... to ... do
    ...
od od;

> q:= eval(max(abs(z-2)-2, 1-abs(z-1)),z=x+I*y);
  plots[implicitplot](q, x=0..4, y=-2..2, filled=true, grid=[50,50],
    gridrefine=3, crossingrefine=3, labels=[Re(z),Im(z)], coloring=[red,white]);

Why don't you use Maple rather than Mathematica syntax for your questions?

1) They are not the same, as far as I can tell.

> [FromMma("Arg[I*Sinh[1 + I]]"),FromMma("Pi - ArcTan[Cot[Abs[Tanh[1 - I]]]]")];

[argument(sinh(1+I)*I), Pi-arctan(cot(abs(tanh(1-I))))]

> evalf(%, 20);

[2.6867724202798433337, 2.6882663475009669513]

2) 
> -I*(-1+z)/(1+z);
  expand(numer(%))/denom(%);

3)  I'm not quite sure what you want, but perhaps

> with(plots):
  conformal(exp(z),z=-5 .. ln(2)+2*Pi*I, axes=box, numxy=[10,60]); 

I just copied the folder C:\Program Files\Maple 12\toolbox\GlobalOptimization into the corresponding Maple 13 folder, included ""C:\\Program Files\\Maple 13\\toolbox\\GlobalOptimization\\lib" in libname, and it seems to work fine.

1) I don't know what you mean by "distribute the result globally" or "globalize".  fsolve produces a result: it's up to you to assign the result to a variable, if that's what you want.  There's certainly no need to cut and paste.

2)  I really don't think fsolve does that.  Any variables that come out of fsolve are the same as the ones that go in.
 

Can you give us examples?

sort(L, `>`);

Rather than use dsolve, you might be able to use odetest and solve(identity...).  For example: find solutions of
  `y''`(x)= y(x)^2/x
of the form
y(x) = a*x^p

> de := diff(y(x),x$2) = y(x)^2/x;
  form := y(x) = a*x^p;
  R:= odetest(form, de);
  
R := a*x^(p-2)*p^2-a*x^(p-2)*p-a^2*x^(2*p-1)

This has to be identically 0 in order to have a solution.

> solve(identity(R,x));

{a = 0, p = p}, {a = 2, p = -1}

So the two solutions of that form are y(x) = 0 and y(x) = 2/x.

Note that ln(cosh(t)) = t + ln((1+exp(-2*t))/2).  Expressing it this way avoids the overflow (there's still an underflow, but it's harmless here: exp(-2*t) is taken as 0 when t is large enough).

So

> zx := 2000 - (14.9165/9.81)*(60*9.81*x/14.981 + ln((1+exp(-2*60*9.81*x/14.981))/2));
 

There is a LieAlgebras package that looks like it may have something like what you want.  See the help page ?DifferentialGeometry,LieAlgebras.  In particular that page has a link to tutorials on using this package.

It is unlikely that solve could make much headway on this problem: it is for closed-form solutions, and that is not likely to be available for expressions involving Bessel functions.  Actually, since your expression involves floats, I think solve passes the problem to fsolve, which looks for numerical solutions.  Unfortunately, due to the huge numbers involved, your expression is very bad from a numerical point of view.  It can't even be evaluated for omega = 0, for example, because this would involve the Bessel functions at huge imaginary numbers.  However, I was able to find some approximate solutions.  For example:

> Digits:= 50:
  fsolve(fm,omega = 0.1162e17 .. 0.1163e17);
                                                             17
      0.11620006344909625656144701396794086246099468957363 10

> fsolve(fm,omega = 0.11621e17 .. 0.11622e17);
                                                             17
      0.11621499999998203909780473614553381926280573214772 10

> fsolve(fm,omega = 0.11622e17 .. 0.11623e17);
                                                             17
      0.11622499999996876934361503106363551591883445176429 10

Actually it seems there is a huge number of solutions in this region.

> plot(1e-30*eval(fm,omega=0.11620006e17+t),t=0..10000);

You have an irreducible polynomial of degree 6 in z.  It has no solution in radicals: the Galois group is S(6).  The RootOf stands for a root of the given polynomial (which is just 1/3 times your polynomial, with z replaced by _Z.   This is the best that can be done.  Of course, for numerical values of c, you could get numerical approximations of the solution using fsolve.  Or you could try for series solutions.  For example, given that z=1 is a solution when c=-20,

> z1:= 819*z^6+(-396+396*c)*z^4+306*z^3*c+(63*c^2-126*c+63)*z^2
    +(90*c^2-90*c)*z+6*c^3-9*c^2+18*c-6;
  S:= solve(z1,z);
  series(S, c=-20);

series(1-1/24*(c+20)-1/9216*(c+20)^2-1/196608*(c+20)^3-7/452984832*(c+20)^4-9397/7044820107264*(c+20)^5+O((c+20)^6),c=-20,6)

Or you could plot the rather interesting-looking curve z1=0.

> algcurves[plot_real_curve](P,c,z,axes=box,labels=[c,z]);

You need to give the name of the variable as well as the range.  Thus:

> plot(tst, p=Pi..2*Pi);

> B:= LinearAlgebra:-Copy(A);
  for i from 1 to 5 do B[i,i]:= 1 end do:

If these are column vectors with m entries,

> W := Matrix(m, h, (i,j) -> w[j][i]);