Here is one solution.  I'm not confident that this will work on Maple 12.01, which I don't have immediately at hand. I'm using a development version and the event mechanism has changed a bit, so if this doesn't work, let me know. Because there isn't a way, I believe, to assign to an Array in the integrator event cycle, I force a halt each time the condition is met, save the results, reinitialize, and continue integrating until reaching the end.

deq := { diff(y(t),t,t)+sin(t)*diff(y(t),t)+y(t)^2-5*y(t)=0 }:
ic := { D(y)(0)=0, y(0)=1 }:
dsys := deq union ic:
integ := dsolve(dsys
                , numeric
                , events=[NULL
                          , [[y(t)-4=0, diff(y(t),t) > 0], halt]
                         ]
               );
Tf := 10:
_Env_in_maplet := true: # trick to quiet the warnings

data := table():
cnt := 0;
do
    result := map(rhs,integ(Tf));
    if integ('eventstop') <> 0 then
        # The halt was caused by the event (not by reaching Tf), so 
        # save the returned value in the data table.
        cnt := cnt+1;
        data[cnt] := result;
    end if;
    # The independent variable (time) is always first in result.
    # Exit loop if we've reached the end, otherwise reinitialize, using
    # result.  Without reinitializing, the integrator halts immediately
    # because the halt condition is fulfilled.
    if Tf <= result[1] then break; end if;
    integ('initial' = result);
end do:

# Convert table to rtable.
data := rtable(1..cnt, data);

f := 4/Pi*sum((sin(2*n-1)*x)/(2*n-1),n=1..infinity);
combine(f,ln);
                                        x

You might try issuing the Maple command

interface(labeling=false);

I'm not sure how to do that in the Maple toolbox for Matlab.

f := (x,y) -> `if`(x=0 and y=0, 0, x*y/(x^2+y^2));
 

While assign does what you asked, as you progress with Maple you may find that using it is not always the best.  The problem is that by assigning values to the symbols x, y, and z, you can no longer easily access the original equations.   For that reason I generally use eval to to assign an element of a solution to is own symbol.  So, given

sol := {x=1, y=2}:

I would do

xn := eval(x, sol);
                       xn := 1

You might prequalify your numbers to speed things up.  A palindrome with an even number of digits is congruent to 0 mod 11.  Testing for that is very fast compared to the fixed cost of converting to a string, so if numbers in the range you are interested have an even number of digits, you can quickly eliminate 91% of them.

Note that list is not a Maple command; list(30) just returns an unevaluated function.  After assigning that to l, the term l[3], for example, causes Maple to generate a tableref (an implicit table reference).  When you then assign a value to that, Maple allocates a hash table.  You can see this by doing

d := dummy(3):
dismantle(d);
FUNCTION(3)
   NAME(5): dummy
   EXPSEQ(2)
      INTPOS(2): 3
d[3] := Pi:
dismantle(eval(d));
TABLE(4)
   EXPSEQ(1)
   NAME(5): false #[protected]
   HASHTAB(257)
      ......
      HASH(7)
         INTPOS(2): 3
         INTPOS(2): 3
         ....
      ..........................................................................................................
      ..........................................................................................................
      .....................................

It is generally not a good idea to assign directly to elements of a list using index notation.  While that is possible, it is slow and also limited; it won't work for a list with more than 100 elements.

You might specify what you expect.   Do you mean triangulate a polygon?  For example, given the convex quadrilateral ABCD would you want to get ABC and CDA?

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