In general it's a bad idea to use the same name for more than one purpose.  However, you should be able to do something like this:

> m := 5;
   with(Units[Standard]):
   L := 6 * Unit('m');

Both integrations in this example are "non-elementary".  Perhaps the Risch method could be extended to deal with antiderivatives that involve the erf function, but it seems this has not been done in Maple, so we're in the realm of "ad-hoc" solutions rather than complete algorithms.  The procedure `int/rpexp` has code for antiderivatives involving a polynomial times exp of a quadratic, but not exp of higher powers.

No, there is no limit on the number of symbols.  In most (but not all) cases where Maple returns an integral unevaluated, the reason is that there is no way to express the integral in "closed form".  The fact is that, contrary to your experience in a typical calculus course, most integrals can not be expressed in "closed form", and the more complicated the integral the more likely it is to be in this category.

You can give fsolve intervals for the variables.  Actually, the long time seems to be taken by simplify in producing LF[1] to LF[4].  It might be helpful to have some upper as well as lower bounds on the variables, to narrow down the search.  Assuming you want to maximize l, I might also try Maximize in the Optimization package.  But some exploration seems to indicate that this problem is badly behaved, because the arguments of ln become negative (and if you restrict to the region where these are positive, the supremum seems to occur in the limit as one of these arguments approaches 0)

I'm not sure I understand your terminology.  There may be a language barrier here.  I would say the three representations of a surface are

1. implicit, e.g. f(x,y,z)=0
2. explicit, e.g. z = Z(x,y)
3. parametric, e.g. x = X(s,t), y = Y(s,t), z = Z(s,t)

The explicit form could be considered as a special case of the implicit (with f(x,y,z) = z - Z(x,y)) or as a special case of the parametric (with s = x and t = y). 

In general it is not a simple problem to convert from one form to another.   For the case of curves rather than surfaces you could use parametrization in the algcurves package to convert from an implicit to a parametric representation (under certain conditions).   I don't know of an equivalent for surfaces.  Of course, for simple examples such as yours a human can produce a parametric representation using geometric insight: the surface is a sphere, and parametric representations of a sphere are well-known.  Or somewhat more generally, if you can guess the general form of a parametric representation you can try solving for the coefficients that would make it satisfy the implicit equation.

Somehow the procedure atataaV seems to be producing a Matrix that contains an unevaluated n.  Since you haven't told us the definition of atataaV, I can only speculate on how that occurs, but one way would be using quotes 'n'.  Now assume actually assigns a new value (sometimes denoted n~) to n, so that Equal sees two versions of the Matrix, one with n unevaluated and the other with n~, and concludes that these are different.  I suspect that inserting an eval in an appropriate place in the definition of atataaV would fix this problem.

All you really need to do is to use plots[loglogplot] instead of plot, and start m from some positive value such as 0.1 instead of 0.

While I have your attention, you may find the following hints helpful:

1) with(plots) should come after the restart, not before: restart wipes out anything that has been done before it.

2) It's very inefficient to build up a list by appending one item at a time, as you've done.  It's much better to use seq.  Thus you could have written

> restart; with(plots):
> plotlist:= [seq(loglogplot(sum(m^i*exp(-m)/i!, i = 0 .. c), m = 0.1 .. 40), c = 0 .. 40)]:
> display(plotlist);

Your k can be simplified:

k:=diff((q*l^(1/3)*q)*(1-l)^(1/2)+(1-a)*(1-q*l^(1/3))*q*(1-l)^(1/2)+a*(1-q*(1-l)^(1/2))*q*l^(1/3), l);
normal(k);

                    (2/3)      (2/3)                1/2
             q (-3 l      + 3 l      a + 2 a (1 - l)   )
             -------------------------------------------
                             2/3        1/2
                          6 l    (1 - l)

Neglecting the case q=0, you want those values of l that make -3*l^(2/3)+3*l^(2/3)*a+2*a*(1-l)^(1/2) = 0.  Note that this does not involve q.

There seems to be a problem with your second equation and fourth boundary condition.  
Note that this equation is singular at r=0, so specifying boundary conditions there may be tricky.  

>   de2:= diff(r*diff(phi(r),r),r) - m^2 *(r*phi(r)-n/m^2)=0;
     S:= dsolve(de2,phi(r));

                           / 2        \
           /d        \     |d         |    2 /            n  \
    de2 := |-- phi(r)| + r |--- phi(r)| - m  |r phi(r) - ----| = 0
           \dr       /     |  2       |      |             2 |
                           \dr        /      \            m  /


  S := phi(r) = BesselI(0, m r) _C2 + BesselK(0, m r) _C1

               n Pi StruveL(0, m r)
         - 1/2 --------------------
                        m

series(rhs(S), r, 2);
             /      /    m r         \    \            2
             |_C2 + |-ln(---) - gamma| _C1| - n r + O(r )
             \      \     2          /    /


In order for the solution to be finite at r=0 you need the constant _C1
 to be 0.  But then the derivative at r=0 is -n, not 0.

The first problem is that the plane command needs points (as defined in geom3d), not Vectors.  So you might try

> with(geom3d):
   point(A, [a,b,c]):
   point(B, [e,f,g]):
   point(C, [h,i,j]):
   plane(p1,[A,B,C]);

Unfortunately this doesn't work: I get an error message

Error, (in geom3d:-plane) unable to determine if A, B and C are distinct

The problem is that geom3d doesn't subscribe to the usual Maple philosophy that things are treated as
distinct unless they are known to be equal. In principle it sounds good: you can't make a plane with
three points unless those points are distinct (and not just distinct but not collinear). But this makes it
very difficult to work with "generic" points.

If you want the equation of a plane through three points, and know a little vector algebra, I think it's better to
use the VectorCalculus package.

> restart; with(VectorCalculus):
   A:= <a,b,c>: B:= <e,f,g>: C:= <h,i,j>:
   N:= (B-A) &x (C-A):  # this is the normal to the plane
   p1:= N . <x,y,z> = N . A;

If the curves are given parametrically, one way to join them by a surface is to interpolate linearly.  For example:

Curve1:= [t + t^2, 1 - t^2, t, t = -1 .. 1];
Curve2:= [t - t^2, 1+t^2, 2, t = -1 .. 1];
with(plots):
Surf:= [op(expand(s*Curve1[1..3] + (1-s)*Curve2[1..3]))];
display([spacecurve(Curve1, colour=red, thickness=3), 
  spacecurve(Curve2, colour=blue, thickness=3), 
  plot3d(Surf, s=0..1, t=-1..1)]);

Not quite sure what you want to animate, but perhaps something like this?

with(plottools): with(plots):
 nn:= 19: A:= <1.2,0.5; -0.3,1.2>: f:= x -> <cos(x), sin(x)>:
 points:= evalf([seq(f(2*Pi*i/nn),i=1..nn)]):
 targets:= [seq(map2(`.`,A^i,points),i=0..4)]:
 ellipses:= [seq(plot([((A^i).f(x))[1],((A^i).f(x))[2],x=0..2*Pi]),i=0..4)]:
 arrows:=[seq([seq(arrow(targets[i-1][j],targets[i][j],difference=true,colour=HUE(j/nn)),    
    j=1..nn)],i=2..5)]:
 display([seq(display(arrows[i],ellipses[i],ellipses[i+1]),i=1..4)],insequence=true, 
    scaling=constrained); 

> ListTools[Split](`=`,a,class);

  [], [23, 45, 74, 34, 25, 76],

        [55, 67, 78, 34, 65, 12, 23, 45, 64, 23], [12, 13, 15],

        [99, 98, 97, 59]

Just to expand on pagan's answer:  assuming basically does this:

1) replace names in the given expression by versions of them on which assumptions have been made

2) evaluate the expression

3) substitute the original names for the "assumed" versions in the result.

In this case, the result of the RandomVariable function is just a name _R.  The details of the distribution of _R are hidden away somewhere as "side effects" of RandomVariable, and are out of reach of the last step of assuming, so they still contain the assumed version of q.

As a general rule, there's no point in using assuming when it's not needed.  RandomVariable and Mean are quite happy to use the distribution Geometric(q) without knowing that 0 < q <= 1, even though the result doesn't make sense if q is outside this interval.  The place where it might make sense to use assuming is in the solve command (although in this case it's not needed).

You might also try fprintf, which does accept Arrays.  See the help page ?rtable_printf for various options.

For example:

> fprintf("c:/test.txt","%d",a);
   fclose("c:/test.txt");

which produces a file that looks like this:

2 3
4 5