You might find the following plot somewhat nicer looking:

pts := [$1..100]:
DynamicSystems:-DiscretePlot(pts,map(numtheory:-pi,pts),style=stair);

Maple doesn't generate a closed-form expression for W(t).  That is, more or less, the whole point of a numerical integrator.  You can inspect the procedure that generates it, however, that generally is not terribly informative.

dsol := dsolve({diff(x(t),t)=sin(t), x(0)=0}, 'numeric'):
interface(verboseproc=2):
print(dsol);

The idea is simple enough.  As mentioned earlier, the circuit has two pairs of complex conjugate purely imaginary poles, that is clear from simple inspection.   The response, then, has the form

A1*sin(w1*t + phi1) + A2*sin(w2*t + phi2)

Think of each angle (w1*t + phi1, etc) controlling the angle of a hand on a clock. The maximum possible values occurs when both hands are straight up, that value is A1 + A2.  That is guaranteed to occur at some time if w1/w2 is irrational (the hands might not meet precisely, but they will be arbitrarily close).  If w1/w2 is rational then whether both hands meet at the top depends on the ratio and their starting angles (phi1 and phi2).

Now, the trick is to solve for A1 and A2.  In another thread you have had some advice on how to do that.  The real trick is reducing the result to a relatively nice form.  Here is a reasonable presentation:

Let

    v = x + y + z
    w^2 = v^2 - 4*x*y
    x = Lm*Cs
    y = Lp*Cp
    z = Lm*Cp

Then

  Apk = A1+A2 = I0*sqrt(2)/8*sqrt(x/y)/w/Cs*((v+w)*sqrt(v-w) + (v-w)*sqrt(v+w))

With some work you can probably reduce this to a decent expression.  The stationary value vs Cp occurs when

  Cp = Cs*Lm*Lp/(Lm+Lp)^2

at that value, all the energy will transfer to Cs (at some point in the response).

The general form of the characteristic polynomial has two pairs of conjugate, purely imaginary roots.  Let the corresponding natural frequencies be w1 and w2.  The voltage at Cs then has the form

  Vcs = V1*sin(w1*t+phi1) + V2*sin(w2*t+phi2)

with V1 and V2 positive constants proportional to the initial current and the circuit parameters.  Provided w1/w2 is not a rational number, Vcs_max = V1 + V2.  Alas, the form of Vcs_max is rather complicated.  If you need more help, I'll post a worksheet showing how you can find it.

Note, by the way, that for any values of Lm, Lp and Cs, there is a value of Cp such that all the energy is transferred to Cs, at some time.

That's a rather general question.  What sort of exploration did you have in mind?

Because there are no dissipative elements (resistors), there is no final equilibrium state.  However, that also means that no energy is lost.  Consequently the maximum value across Cs occurs when it stores all the energy, provided that occurs.  That is a trivial computation. Whether it occurs (all energy transferred to Cs) depends on the ratios of the time-constants (I believe).

I do not believe that they implemented their algorithm in Maple. Rather, they compared their algorithm to a symbolic analysis using Maple. That is, for the Maple side they generated the nodal equations and fed them to Maple's solver.

You can access the paper here www.ee.ucr.edu/~stan/papers/iccad97.pdf

For n points, with m collinear I expect the answer to be

binomial(n - m, 2)*m + binomial(m, 2)*(n - m) + binomial(n-m, 3)

An alternative is to use map, which returns a list (because its second argument is a list):

map(op, [1,3,4], a);

                             [2, "john", sin(57)]

Here's a clever approach.  It assigns a package that reassigns the one Maple postfix operator (!).

tailevalf := module()
option package;
export `!`;
    `!` := proc()
        evalf(_passed);
    end proc;
end module:

> with(tailevalf);
                                      [!]

> sin(Pi/3)!;
                                 0.8660254040

> unwith(tailevalf):
> 5!;
                                      120

This forum is for Maple users; I'm guessing you aren't using Maple.  Get all x terms on one side, factor out the x, then divide by the coefficient.  For example:

 x/2 + x/6 = 1

 x*(1/2 + 1/6) = 1  # factor out the x

 x*(2/3) = 1           # combine fractions

 x = 3/2                 # divide by coefficient of x

See the help page ?worksheet,documenting,styles.  The first subsection explains how to save a style as the default.

Did you try collect(T1,x).  Is that what you want?

You might scale to make a polynomial and then use PolynomialTools:-CoefficientList:

PolynomialTools:-CoefficientList(expand(T1/x^(min(0,ldegree(T1,x)))),x);
[2/3 q[1], 2 + 2/3 q[2], Pi I + 2 q[1], -p[0] + Pi q[1] I + 2 q[2],

    -p[1] + Pi q[2] I, -p[2], -p[3]]

I doubt there is such an expansion.  It isn't too hard to see why.