@tomleslie That's very nice.  I did not think of applying parse() in this context.  Your construction can be handy for displaying the value of one variable in a title.  I don't quite see how to extend it to more than one variable but with Carl's answer, that's no longer a pressing issue.

Thank you, Carl, for showing me how.  Before I posted, I had tried all sorts of combinations with typeset() but somehow had failed to hit upon your construction.

Actually the somewhat complex animation code that you have presented obscures the main idea behind the formatting of the title for a single frame which is quite simple. Here it is, should others be interested:

plot(0,
  title=typeset(
    theta[1], sprintf(" = %5.3f   ", 1.000),
    theta[2], sprintf(" = %5.3f   ", 0.123),
    theta[3], sprintf(" = %5.3f   ", 3.333)
  )
);

Addendum: The theta[1], theta[2], theta[3] here may be replaced with theta__1, theta__2, theta__3 with almost equal effect.  See  Acer's comment regarding the difference.

@Preben Alsholm Thansk.  That's good!  In fact the "cheating" is quite mild if we have a rough idea of the shape of the eigenfunction.  The method does not seem to depend very strongly on the guess for omega.

Let's say we are looking for an eigenfunction that is zero at the interval's midpoint (as well as the end points, as they are specified in the boundary conditions).  Then for the approxinate shape we take a cubic:

y(x)=x*(Pi-x)*(Pi/2-x)

Better yet, since we are looking for a solution with D(y)(0) = 1, we multiply by a factor to make it so:

y(x)=2/Pi^2*x*(Pi-x)*(Pi/2-x)

We expect the solution to be omega2 = 4.  We try guesses as wide-ranging as omega2 = 0.1 or omega2 = 100, and in all cases dsolve() finds the correct solution omega2 = 4.

restart;
de := diff(y(x),x,x) + omega2*y(x) = 0;
bc := y(0)=0, D(y)(0)=1, y(Pi)=0;
dsol := dsolve({de,bc}, numeric, approxsoln=[y(x)=2/Pi^2*x*(Pi-x)*(Pi/2-x), omega2=0.1]);
plots:-odeplot(dsol, scaling=constrained);
dsol(0);
dsol := dsolve({de,bc}, numeric, approxsoln=[y(x)=2/Pi^2*x*(Pi-x)*(Pi/2-x), omega2=100]);
plots:-odeplot(dsol, scaling=constrained);
dsol(0);

@Preben Alsholm Ah, that's good.  Now what is the trick to get the other eigenvalues?  Let's say I want to get the eigenvalue omega=2 or equivalently omega2 = 4.  Do you put in an initial guess somewhere?

Hi Preben, I attempted to apply that same method to another eigenvalue problem but it failed:

restart;
de := diff(y(x),x,x) + omega^2 * y(x) = 0;
bc := y(0)=0, D(y)(0)=1, y(Pi)=0;
dsolve({de,bc}, numeric);
    Error, (in dsolve/numeric/bvp) matrix is singular

Is there a way of rescuing the situation?

@Kitonum I had gone as far as op([1,2,1,3,1], ...) to isolate the denominator, but I did not know about applyop() so I couldn't go any farther.  Thanks for offering this solution.  I am happy to have learned about applyop().

@ecterrab Thank you for that extremely compact solution.  I had never thought of using simplify(expression, size) for such a  purpose.

@Axel Vogt Ah, I see that I have misundersood you.  I had wondered a bit why you'd ask that question.

@Markiyan Hirnyk Yes, the winding number is defined with respect to a given point.  I have taken that point to be the origin.  Doing it relative to an arbitrary point is a matter of translating the origin.  I will leave that as a simple exercise to the interested reader.

@Carl Love 

@Axel Vogt There is no need for integration when the curve is piecewise linear.  All we need is to add the angles subtended by the line segments that form the curve.  So we introduce a proc, called subtang(u,v) here, that returns the signed angle at the origin subtended by the line segment connecting the points u, v. If the ray from the origin to a point moving along the line segment from u to v turns counterclockwise, the angle is positive, else it is negative.

subtang := proc(u,v)
  local
    ulen := sqrt(u[1]^2 + u[2]^2),
    vlen := sqrt(v[1]^2 + v[2]^2),
       d := u[1]*v[1] + u[2]*v[2],    # dot product
       c := u[1]*v[2] - u[2]*v[1];    # cross product
  return signum(c)*arccos(d/(ulen*vlen));
end proc:

Now, let P be a list of points in the plane, where a point is a list of two numbers, as in the previous posts in this thread.  We connect the points with line segments, in the order given, to obtain a directed planar piecewise linear curve.  We close the curve by connecting the last point to the first point.  To find the curve's winding number about the origin, we apply subtang() to each of the line segments, add, and finally divide the sum by 2*Pi.  This leads to:

winding_number := proc(P)
  local Q := [P[], P[1]];
  return add(subtang(Q[i], Q[i+1]), i=1..nops(P)) / (2*Pi);
end proc:

Here is a test:

P := [[ -1,-1], [1,-1], [1,1], [-2,1], [-2,-2], [2,-2], [2,2], [-1,2]];

winding_number(P);
                                                   2

plots:-polygonplot(P, style=line, color=red);

@Markiyan Hirnyk Your tests indicate that the timings and memory usage are of the same order of magnitude in all cases.  I conclude that the choice of the method should depend on what best expresses the intent.

@Carl Love That's a tour de force of Maple programming! I find it particularly intriguing that the space dimension does not appear in it explicitly.

@Carl Love Regaring the need for eval() in your comment above, I did just that in the worksheet that I posted earlier today at:

http://www.mapleprimes.com/questions/205350-Natural-Parametrization-Of-Curve

Kitonum and Carl, thank you very much for your very helpful answers to my questions.  Having seen your answers reminded me that I had been told about seq('a[i], b[i]', i= 1..3);  which I had forgotten.

It turns out that there was an extra twist in my application.  Here is what I really wanted:

restart;

c := a[i];

seq('c, b[i]', i=1..3);

which is fine.  This may also be obtained through your other suggestions:

seq(op([c, b[i]]), i=1..3);

or

f := (c,b) -> (c,b);
seq(f(c,b[i]), i=1..3);

The fourth suggestion, however, has different semantics; c is not evaluated:

''c, b[i]'' $i=1..3;

 That's OK, I can pick one of the other three solutions.