@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.

@asa12 It has been pointed out more than once that your original equations are nonautonomous, and therefore a vector field analysis will get you nowhere.

@9009134 The equations in your rewrite.mw are entered in a most disorganized way but I did my best to figure out what you had meant to type, and the way I see it, some things need to be fixed.

1. The conditions
   
   
   are very likely incorrect.  Equations such as yours require boundary conditions, that is, a specification of what the unknowns look like along the domain's edges.  The two conditions above specify phi at the domain's corners.  As such, they are pretty useless.  Perhaps you meant something else.
2. I assume that you mean that the conditions
   
   hold for all -1/2 < z < 1/2.  But if so, then it follows that the conditions shown in the first item above are automatically satisfied, therefore they should be removed.
3. In your attempt to solve the equations you have applied dsolve().  That, however, is good for solving ordinary differential equations.  There is also pdsolve() which is good for solving partial differential equations.  Your equations, however, are neither ordinary nor partial differential equations since they involve the definite integral of the unknown w.  I don't know of a Maple function for solving such equations, so you are really on your own.  If I were doing this, I would try a finite difference scheme.  Writing a program for that shouldn't be too difficult.
4. I am concerned about the wide range of the numerical values of your coefficients.  You have E = 10^11 and f = 10^-12 which span 10^23 orders of magnitude!  Normally one gets rid of such wide span of numerical values through normalizing things with the help of dimensional analysis.  I assumed that the purpose of the bars over the variables in your PDF was to indicate that they had been normalized, but perhaps I am wrong.
   
   You can expect those weird coefficients to be obstacles in implementing any numerical scheme for solving your equations.

@asa12 Your questions seem to be a moving target—they keep changing with each iteration.

Here I will answer the question "how to draw arrow on sphere?"

I assume that you want to draw an arrow which is tangent to the sphere, or, based on your previous questions, you may want to plot a vector field which is tangent to the sphere, like this:

If that's what you want, then have a look at mw.mw

@rit It seems to me that you are mixing up two different concepts:

1. Vector fields
   
   A vector field is something like [f(a,b,c), g(a,b,c), h(a,b,c)].  In Maple you apply fieldplot3d() to plot that vector field as a bunch of arrows, as you can see in the web page on fieldplot3d which you have quoted.  This is independent of any differential equations.
   
   Note that there is no "t" in here.
   
   
2. Differential equations
   
   Associated with that vector field there is a system of differential equations:
   
   da(t)/dt = f(a(t),b(t),c(t));
   db(t)/dt = g(a(t),b(t),c(t));
   dc(t)/dt = h(a(t),b(t),c(t));
   
   In Maple, you apply DEplot3d() to plot the solutions of that system of ODEs.  You will get a bunch of curves.

For instance, the vector field may be [a/r, b/r, c/r], where r=sqrt(a^2 + b^2 +c^2).  You may plot this vector field by applying fieldplot3d().

The associated system of differential equaitons is
          de1 := da(t)/dt = a(t)/r(t);
          de2 := db(t)/dt = b(t)/r(t);
          de3 := dc(t)/dt = c(t)/r(t);
You may plot the solutions curves with DEplot3d().

Now, what else do you want to do?

@Kitonum That's good.  Thanks!