I do not give much for a factor of 2. For real speed one may want to have specific SW for specific tasks. But the price may be high (setting it up, learning). So I would ask: for what? I am satisfied, that Mpl covers that, if I use it as "working horse". It would be more interesting to compare speed beyond double precision (kicking out Matlab). And even then: personally I will spend most of time to feed the problem. Else: use specific SW. Just my 0.02 €.

I would prefer to map a function (also see Preben Alsholm below)

  f:= x -> `if`(1=x, 0, x);
  map(f,A);
                     f := x -> `if`(1 = x, 0, x)


                            [0    0    6]
                            [           ]
                            [5    3    2]
                            [           ]
                            [0    0    0]

A mathematical answer:

This is an implicit definition, valid in Math, if you *know* there is
an explicit, unique solution for f(x) = something (i.e. that f has a
[partial] left inverse).

That is not a good idea for a CAS (and perhaps related to formal proofs)
it is even fals for f = sqrt (having 2 solutions).

One way may be: 

  omega:= 'omega':                        # clear variable
  omega^2 = sqrt(w0^2*(1+((z-zf)/z0)^2)); # set up taks
  [solve(%, omega)];                      # *solve* it (if possible)
  
  omega:= %[1];                           # select (!) the very solution


NB: note that this might make no sense in some context if you leave the
Reals and work in other domains, https://en.wikipedia.org/wiki/Modular_square_root

I do not like *.mw, but sometimes I want it for graphics. Then on WIN I just say "open with ...", using the context menu. And save it as such.

Digits:=15;                 # ~ double precision
eq2:= collect(eq2*N, N);    # avoid divisions
eq4:= collect(eq4*N, N);

{eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9}:
sys:= convert(%, rational); # often better than floats
vars:={C, C1, C2, CT1, CT2, H, N, S, T};

st:=time():

solve(sys, vars):          # solve symbolically
evalf[60](%):              # care for numerical errors
evalf(%);                  # show in desired precision

# gives 4 solutions in 10 (?) seconds
time() - st;


Edited: with some more care I get 11 solutions. The first 4 are by roots of 2 
polynomials of degree 2, the remainings are roots of a polynomials of degree 7

sols:={solve(sys, vars)};
nops(sols);
                                  4
sols[1..3];
map(allvalues, %); 
evalf[60](%): 
S13:=evalf(%); #nops(%); 

sols[4]:
{allvalues(%)}:
evalf[60](%): 
S4:=evalf(%);  #nops(%);

S13 union S4; nops(%);

                                  11

The output is somewhat lengthy.

Have not checked whether all are correct solutions (may be Digits = 15 is not enough here)

(chi^2*(xi^2-1)^2+1-Omega^2)*(xi^2-Omega^2)-xi^2 = 0;
eval(%,   chi^2=1/1200);

plots[implicitplot](%, xi=0 .. 2, Omega=0..2);

Remark beyond pretty presentation:

1) Perhaps I would use Modulo(5^`1000`, 29) in classical interface for
typing (but forgot how to resolve the name `1000`)

2) For efficient evaluation I still have in mind Fermat's little theorem

  p:=29;
  1000 mod (p -1); # prepare for little Fermat
  5^% ;            # considerably smaller 
  % mod p;         # now use Fermat

                                  23

Wikipedia, "little Fermat"

Maple has a help and handbooks. You may wish to have a look.

I think you do not want to solve for theta_n, but you actually want the
tan of (and then invert by arctan), because of your range.

eval(EQ, theta_n = arctan(x));
f:= unapply(%, x,  n,   omega_a, v, x_l, C_a, Z_0); # solve for x

Now plotting some examples (using your parameters) indicates that a good
guess is given by the intermediate value at say x=+-10

  initial_x:= (F(10) - F(-10))/2; # F = f filled with data

Which means to use arctan(10) ~ 1.5 (avoiding calculations).

I get a speed up of ~ 4 times, only, but that may indicate a way and
using compiled versions you may gain a factor of 10.

Note that you can set t=1 (by combining it into a).

Since cos(r) is the real part of exp(r*I) you can write your integrand as
Re( cos(xi*x)*exp(-xi^2*c) ) for c = b^2-I*a. Now use the assumptions:

  Int(cos(xi*x)*exp(-xi^2*c),xi = 0 .. infinity);
  eval(%, c = b^2-I*a);
  value(%) assuming 0<a,0<b,0<x; 
  Re(%);
  simplify(%, {b^2-I*a = c});

                   infinity
                  /
                 |                             2
                 |          cos(xi x) exp(-c xi ) dxi
                 |
                /
                  0



                                             2
                            / Pi \1/2       x
                     1/2 Re(|----|    exp(- ---))
                            \ c  /          4 c

x^3-3*x^2-24*x+8:
f:=unapply(%, x);

X:=[-2,4,6, 9];

map(f, X);
                               3      2
                    f := x -> x  - 3 x  - 24 x + 8

                         [36, -72, -28, 278]

I suggest to take any serious book on algebraic Geometry. And rephrase after reading. Sorry for being harsh.

I do not understand why teaching such is obfucated by using a CAS ...

Expand the power and divide by x, giving x^2 + 6*x + 12 + 16/x.
For x towards 0 only the last term matters. Now care for sign,
in case of doubts take some values to understand.

And if "handwaving" is allowed like in school you can do shorter.

s:=proc(n::posint) {seq(j, j=1..n)}; end proc;

s(4); 
                             {1, 2, 3, 4}
whattype(%);
                                 set

Just some remarks aside, looking at some of your questions: it would make
sense to learn some basics Maple 'syntax'.

For example try to avoid f(x):= ... as this does not define a function
(though it is used in Math sometimes). Either f:= x -> (some expression)
or by 'unapply' is a much better way (as you used later).

Also never hesitate to use brackets: even if you might think it is 'clear'
what is meant it may be not (cf Carl's question) and a machine may moan.

For old versions there are guides like 'first steps', they are dated and
may not be fine for the 'standard interface'. But I think they are "the"
right way to learn basics:

http://www.maplesoft.com/view.aspx?SF=94/LearningGuide.pdf