About is, its capabilities and limits, and the relation of CAS with theorem provers, you may see e.g. these threads:

Theory-Prover-Capability-In-Maple

Is--If--With-Assuming

Type-Complexcons

May be that it sounds a bit confusing as it was not answered yet, but it is quite simple if you want to use Maple:

W := Q*qA - P*(qA - (qA-qE)/2) - (1/2)*k1*(qA - (qA-qE)/3)^2 - (1/2)*k2*(qA - 2(qA-qE)/3)^2:
dWqA := diff(W,qA);
dWqE := diff(W,qE);
                                /2 qA    qE \
                           2 k1 |---- + ----|
                                \ 3      3  /
         dWqA := Q - P/2 - ------------------ - k2 (qA - 2/3)
                                   3


                                      /2 qA    qE \
                                   k1 |---- + ----|
                                      \ 3      3  /
                   dWqE := - P/2 - ----------------
                                          3

solve({dWqA,dWqE},{qA,qE});

        6 Q + 3 P + 4 k2         12 k1 Q + 6 k1 P + 8 k1 k2 + 27 k2 P
  {qA = ----------------, qE = - ------------------------------------
              6 k2                             6 k1 k2

        }

Actually, shipped Maple includes diverse tensor facilities including the older tensor package (cf. ?tensor) and the newer Tensor package (cf. ?DifferentialGeometry,Tensor). Besides, there are several user contributed packages. Each of them has a different approach to the subject, and which is best (if any) depends on the details of the intended computation.

A way to get that output form is using evaln:

n:=3:
a:=(r)->q[r](t):
x[1]:=Vector[row](n,a); 
w[1]:=x[1];                  
for j from 2 to 4 do
w[j]:= Vector[row]([evaln(w[j-1]),x[j]]);  
end do;

                 x[1] := [q[1](t), q[2](t), q[3](t)]

                 w[1] := [q[1](t), q[2](t), q[3](t)]

                         w[2] := [w[1], x[2]]

                         w[3] := [w[2], x[3]]

                         w[4] := [w[3], x[4]]

Perhaps, it might be simpler to try another approach, go back to the ode for which the starting equation is a solution, inserting values for the parameters, and solve numerically for the IC:

eq:=v*t=(S0-St)+Km*ln(S0/St):
PDEtools:-dpolyform(subs(St=St(t),Km=11, S0=1,v=1,eq), no_Fn);
   d              St(t)
  [-- St(t) = - ----------] &where
   dt           St(t) + 11

        [-St(t) - 11 <> 0, t - 1 + St(t) <> 0]

ode:=op([1,1],%);

                           d              St(t)
                    ode := -- St(t) = - ----------
                           dt           St(t) + 11

sn:=dsolve({ode,St(0)=1},numeric):
sn(1);
                 [t = 1., St(t) = 0.919783724972999]

Or write a series solution:

 
dsolve({ode,St(0)=1},St(t),series);

                        11   2     11    3      143     4
  St(t) = 1 - 1/12 t + ---- t  - ------ t  + --------- t  +
                       3456      165888      286654464

            2629      5      6
        ------------ t  + O(t )
        206391214080

If you want a series for small values of the LambertW argument, u say, you can do something like:

eq:=v*t=(S0-St)+Km*ln(S0/St); 
                                               S0
                  eq := v t = S0 - St + Km ln(----)
                                               St

s:=solve(eq,St);
                                      -v t + S0
                               S0 exp(---------)
                                         Km
                 s := LambertW(-----------------) Km
                                      Km

s1:=subs((op([1,1],indets(s,specfunc(anything,LambertW))))=u,s);
                         s1 := LambertW(u) Km

series(s1,u=0,3);
                                    2      3
                         Km u - Km u  + O(u )

This has changed in Maple 14:

convert(LambertW(x),Int):
series(%,x):
Error, (in series/int) unable to compute series

about(x);
x:
  nothing known about this object

For me, Maple 14 Standard

eqn := expand((x+y)^2000) = 0;

prints all the terms and there is no apparent difference with or without

interface(elisiontermsthreshold = infinity);

Maple 12.02 and Maple 13.02 Standard also prints all the terms, but much slower without the elision statement. On the other hand, I get a warning message here:

eqn := expand((x+y)^10000) = 0;

                 [Length of output exceeds limit of 1000000]

with or without the elision statement.

It may be worth noting that Maple 14 does not provide the "traditional" Classic GUI Help Browser and Classic help pages, but only the Standard GUI Java Help Browser (to my regret).

I understand that the XML help files are filtered for reading them as ascii text within the console UI, when  interface option helpbrowser is set to text. Otherwise, when set to  standard, the Java browser is called.

You can plot with:

plot3d(k,  -1 .. 1,  -1 .. 1);

Maple 14 includes "experimental" routines for exact solutions for PDEs subject to boundary conditions, see ?updates,Maple14,de. They provide mode expansions for some problems, but not this one. Trying it, just gives NULL. This is a trace of the methods it tries before that:

PDE:=diff(u(x,t),t$2)+diff(u(x,t),t)-diff(u(x,t),x$2)=sin(Pi*x/l):
BC:=u(0,t)=0,u(l,t)=0:
IC:=u(x,0)=0,D[2](u)(x,0)=0:

infolevel[pdsolve] := 2:
pdsolve([PDE,BC,IC]);
* trying methods for class "_Fn" for 2nd order PDEs
   -> trying "linear_in_xt"
   -> trying "BC_equal_0"
* trying methods for class "_Cn_cn" for 2nd order PDEs
* trying methods for class "Wave" for 2nd order PDEs
   -> trying "Cauchy"
   -> trying "SemiInfiniteDomain"
   -> trying "WithSourceTerm"
* trying methods for class "Heat" for 2nd order PDEs
   -> trying "SemiInfiniteDomain"
   -> trying "WithSourceTerm"
* trying methods for class "Laplace" for 2nd order PDEs
   -> trying a Laplace transformation
* trying methods for class "Fourier" for 2nd order PDEs
   -> trying a fourier transformation
* trying methods for class "Series" for 2nd order PDEs
   -> trying "sincos"
   -> trying "2"
   -> trying "3"
   -> trying "4"
   -> trying "WithSourceTerm"
* trying methods for class "Generic" for 2nd order PDEs
   -> trying a solution in terms of arbitrary constants and functions to be adjusted to the given initial conditions

You could enter the speed as:

v:=3*10^6*Unit('m')/Unit(s);
                                3000000
                           v := ------- [m]
                                  [s]

You can do this:

eqn:=diff(x(t),t)+diff(y(t),t,t)+x(t)=0:
PDEtools:-Library:-DifforderJet(PDEtools:-ToJet(%,[x,y]),[x[t],y[t]]);

                            [x = 1, y = 2]

For instance:

ode:=diff(y(x),x$3)-2*diff(y(x),x$2)+x^2*diff(y(x),x)+sin(x)*y(x)-1=0:
with(DEtools):
convertsys(ode,{}, y(x), x, Y);

  [[YP[1] = Y[2], YP[2] = Y[3],

                              2
        YP[3] = 1 + 2 Y[3] - x  Y[2] - sin(x) Y[1]],

                                              2
                             d               d
        [Y[1] = y(x), Y[2] = -- y(x), Y[3] = --- y(x)], undefined, []
                             dx                2
                                             dx

        ]

You can program the transformation:

f:=proc(x) local a,b,c,d;
a:=op(1,x);
b:=op([2,1,1],x);
c:=op([2,1,2],x);
d:=op([2,2],x);
(a^(1/d)*b + a^(1/d)*c)^d;
end proc:

f(10.17594299*(-.2450501048*x+.2692858295)^2.824978580);
                                                 2.824978580
                 (-0.5570846540 x + 0.6121809386)