Did you try it with the compiler option selected? 40 seconds seems slow.  On my 920 it runs in less than half a second.

You haven't assigned accel.  The call to UpdateAllVelocities raises an error because the passed arguments don't match the declaration.  Assign accel a Matrix and it works.

The usual way to do this is to call evalf.  Better yet, express the sum in an inert form (using ?Sum) and then call evalf. That is

 S := Sum(ln(x)/((x)^(2)+3*x+2),x=1..infinity):
 evalf(S);

Alas, Maple is not able to numerically evalute the sum, at least not without a little help. Try the following

S1 := applyop(convert,[1,2],S,parfrac,x);
                           infinity
                            -----
                             \            /  1       1  \
                     S1 :=    )     ln(x) |----- - -----|
                             /            \x + 1   x + 2/
                            -----
                            x = 1

S2 := applyop(expand,1,S1);
                              infinity
                               -----
                                \      /ln(x)   ln(x)\
                        S2 :=    )     |----- - -----|
                                /      \x + 1   x + 2/
                               -----
                               x = 1

S3 := map(op(0,S2), op(S2));
                       /infinity      \   /infinity          \
                       | -----        |   | -----            |
                       |  \      ln(x)|   |  \      /  ln(x)\|
                 S3 := |   )     -----| + |   )     |- -----||
                       |  /      x + 1|   |  /      \  x + 2/|
                       | -----        |   | -----            |
                       \ x = 1        /   \ x = 1            /

At this point it is easier to do a bit of manual work.

S3 = (ln(1)/2 + ln(2)/3 + ln(3)/4 + .... ) - ( ln(1)/3 + ln(2)/4 + ln(3)/5 + ... )
   = ln(1)/2  + (ln(2)-ln(1))/3 + (ln(3)-ln(2))/4 + ... 
   = sum(ln(1+1/k)/(k+2), k = 1..infinity)

so

S4 := Sum(ln(1+1/k)/(k+2), k = 1..infinity):

evalf(S4);
                          0.6017824584

As a check

evalf(add(op(1,S),x=1..10000));
                                 0.6007616145

Here's one approach.

Clean := proc( Opts :: list, P :: list)
local opts, k, nm;
    ASSERT(nops(Opts) = nops(P));
    opts := Opts;
    for k to nops(P) do
        if P[k] = -1 then
            nm := opts[k];
            opts := map(L -> `if`(L::list
                                  , remove(`=`,L,nm)
                                  , L
                                 )
                        , opts);
        end if;
    end do;
    return opts;
end proc:

Clean([Joe, [Mike, Abraham, Leon, Joe], [Terresa, Cody, Joe]], [-1,2,2] );
                 [Joe, [Mike, Abraham, Leon], [Terresa, Cody]]

While it isn't impossible, it is difficult and probably not worthwhile to do that.  As assigned, the expression ALL does not contain L1, but rather the value of L1.   To do what you want you could use forward quotes when assigning ALL:

ALL := '[L1,L2]';

Now ALL contains the names L1 and L2.  With some care and use of the uneval parameter modifier (see ?parameter_modifiers) and the ?eval procedure, you can inspect and manipulate both the name and its value.  However, doing so is tricky.

My suggestion is that you don't do this.  A simpler approach is to use a table

T[L1] := [1,2,3]:
T[L2] := [2,3,4]:
T[ALL] := [L1,L2]:

Then you could do

for L in T[ALL] do
   if StringTools:-RegMatch["1", L] then
        T[L] := map(`+`, T[L], 1);
   end if;
end do;

Here's an approach using ?rsolve, Maple's recurrence equation solver

G:=s->q*(1-p*s)^(-1);
                                            q
                               G := s -> -------
                                         1 - p s


req := g(k+1) = G(g(k));
                                               q
                         req := g(k + 1) = ----------
                                           1 - p g(k)


rpoly := simplify(denom(rhs(req))*req);
                     rpoly := (-1 + p g(k)) g(k + 1) = -q


rsol := simplify(rsolve({rpoly, g(0)=s},g(n)));
                         1/2      n               1/2             n
rsol := (s %2 (1 - 4 q p)    + s 2  %1 (1 - 4 q p)    + s %2 - s 2  %1

                 (1 + n)         /                      1/2               n
     - 2 %2 q + 2        %1 q)  /  (-%2 + %2 (1 - 4 q p)    + 2 p %2 s + 2  %1
                               /

        n               1/2    (1 + n)
     + 2  %1 (1 - 4 q p)    - 2        p %1 s)

      /         q         \n
%1 := |-------------------|
      |                1/2|
      \-1 + (1 - 4 q p)   /

      /         2 q        \n
%2 := |- ------------------|
      |                 1/2|
      \  1 + (1 - 4 q p)   /


simplify(eval(rsol,n=0));
                                       s

simplify(eval(rsol,n=1));
                                        q
                                   - -------
                                     s p - 1

simplify(eval(rsol,n=2));
                                 q (s p - 1)
                                --------------
                                -1 + q p + s p

Look at the value of SEqn, it has the form

{ (1+x1)*exp(x1) = Int( f(v), v = k*Eg1 .. infinity), ... }

where k is a number.  fsolve won't handle those.  It returns unevaluated.

Here's the basic idea

plots:-animate(plots:-pointplot, [[i,i^2]], i=[seq(1..100)], symbol=circle);

For data generated in a loop the following should work

n := 100:
for i to n do
    data[i] := i^2;
end do:

plots:-animate(plots:-pointplot, [[k,data[k]]], k=[seq(1..n)], symbol=circle);
Error, points cannot be converted to floating-point values

Alas, there is a bug in `plots/animate`, it isn't evaluating the result after plugging in values. I'll submit an SCR against that.

Apparently one of the arguments must be integer[8] to work properly:

for i in [1,2,4,8] do
    for j in [1,2,4,8] do
        A := Array([1], 'datatype'=integer[i]);
        B := Array([1], 'datatype'=integer[j]);
        C := A +~ B;
        if C[1]::integer then
            printf("%d %d\n", i, j);
        end if;
    end do;
end do;
1 8
2 8
4 8
8 1
8 2
8 4
8 8

I suspect this was done to reduce overflow issues. The safer thing might be to create the output as datatype=anything, but that isn't ideal.

You might also consider using ?iquo, however, that doesn't work with non integer numeric arguments.

One approach is to use the ?difforms package (a rather old package, but still useful for some tasks).

y := r*sin(theta);
                               y := r sin(theta)

with(difforms):
defform(r=scalar,theta=scalar);
d(y);
                    sin(theta) d(r) + r cos(theta) d(theta)

Fleshing out Alec's suggestion:

evalf[30](eval({seq(Eqn(j),j=1..n)}, EV));

It's not just more elegant, it is O(n) instead of O(n^2).

How are you computing the sequence?  Here's one approach to plotting this

x_next := proc() global x; x := 9*sqrt(x+6) end proc:
x := 2.30:
plot([[1,x],seq([k,x_next()], k=2..30)]);

As for the matrix problem, the issue there is that Maple is computing the exact result, so the fractions contain large integers which makes the output look weird.  If you only want floating point values, do

TP := evalf(T)^1000:

Do you want to plot a surface, that is p(t,H) as both t and H vary?  Or do you just want to plot a curve vs t for various H?  For the latter do, for example

deqs := { diff(x(t),t) = h*sin(t) + y(t), diff(y(t),t) = x(t), x(0)=0, y(0)=0}:
integ := dsolve(deqs, numeric, parameters=[h]):
integ(parameters=[1]);
plots:-odeplot(integ, [[t,x(t)],[t,y(t)]], 0..1);

Are you using the ImageTools package that is included with Maple?