I figured out the answer to my own question:

This generalizes nicely to arbitrary non-scalar-valued procedures.  I will explain the details in a Post.

@Carl Love That's illuminating, however I don't understand why you have And instead of Or in the defintion of Nparams.
Consider

f := proc(x, {y:=12})
	x -> x^2;
end proc:

Then Nparams(f) fails as written, but succeeds when we replace And with Or.

@acer that's excellent.  I had tried op(1,q) but that wansn't good.  It had not occurred to me to change q to eval(q).  Thanks for the lesson.

Regarding parameters vs arguments, I thought a bit about it when I posted the question.  To me it seems that the proper usages are "The procedure q has two parameters".  The procedure q takes two arguements" because arguments are passed to q.  But perhaps that's too much of a hair-splitting.

Hello @mmcdara, what you have proposed produces a vector whose first component is D[1](a(u,v)) which is not a very meaningful expression.  We want D[1](a)(u,v) instead.

@one man Getting rid of the denominators in the ODEs is not necessarily a good thing.  It's true that doing that results in neater looking equations, but that does not mean that it improves their solvability.

In the case of the current problem, after removing the denominators, the equations reduce to the deceptively pleasant-looking polynomial ODES:

diff(x(t),t) = z(t)*y(t)^2*(y(t)*z(t)+200)
diff(y(t),t) = -2*x(t)^2*z(t)*(x(t)*z(t)+100)
diff(z(t),t) = 20*(2*x(t)-y(t))*y(t)^2*x(t)^2

ICs := x(0) = 6.066066075, y(0) = -8.537156781, z(0) = 20  # as before

The solution of that system, however, behaves horribly.  The orbit still traces the same curve as before, as it should. However, having removed the denominators, the orbit is no longer a unit-speed curve.  The speed varies very sharply (almost discontinuously) in a couple of places, and as a result, the roots that we are seeking bunch up very close to each other.  So much so that fsolve() has difficulty in finding them.  Here is what the function f3 looks like along the orbit:

The 116 roots are hidden in that mess. Compare that to the corresponding graph in my original post.

By the way, in case you wish to try that yourself, the parameter range for a full orbit is t=0 .. 0.007604467601.

@mmcdara The purpose of my question is pedagogical.  In your solution you refer to L2 and Linfinity norms.  It seems to me that Otttimor does not know what those terms mean.  That''s why I asked him/her for an explanation.

@Otttimor Let's see if this helps.  Consider two lists of numbers:

vals1 := [1,3,5,6,6,6,7,5,3];
vals2 := [1,3,4,5,6,8,7,4,4];

Can you explain how you would want to measure the error/deviation between these lists? Your answer to that will determine the strategy to use in your problem.

By the way, you may plot that data through

dataplot([vals1,vals2], color=[red,blue]);

@Jonas111996 The coefficients of x1(t) and x2(t) on the right-hand sides of your differential equations are things like sin(t*sqrt(-3)) which are explicit functions of t.  That makes your system nonautonomous as Carl has already pointed out, so a phase portrait makes no sense.

Furthermore, due to the presence of the sqrt(-3), those coefficients are complex, making a 2-dimensional phase portrait impossible.

I suspect that you have gone off in the wrong direction in the process of arriving at your differential equations.  If you explain where those equations come from, someone should be able to orient you in the right direction.

You may want to ask the first author whose email is given in the article:

rxyao@sjtu.edu.cn

@acer I agree.  I wouldn't be surprised if a close look at the algorithm reveals better ways of achieving whatever that program is aiming to achieve.

@zenterix Here is what I got when running your file on my Linux machine. 

In the first attempt Maple crashed after running for some half an hour or so.  I montitoried the system's status during the execution and observed that my swap partition, which is a measly 0.8 GB, got filled and apparently that was the cause of the crash.

Then I added a largish (100 GB) swap file and ran the program again.  I didn't time it, but I am guessing that it took about an hour or so to produce the sphere with an arrow.  It turned out that 100 GB of swap space was much more than necessary.  Only a small fraction of it was used.  As I played the animation, however, the memory usage went up until it nearly ate up all of the 64 GB RAM of my computer's memory and the machine started swapping the memory to disk.   At that point I killed the program because the machine had become almost unresponsive.

With your 256GB RAM you may have better luck. 

The lesson learned here is you need a descent sized swap space to run your program.  Type "free" on your Linux command-line to see the size of your swap area.  (The number shown is in kilobytes).  It looks like you will need something like 20 GB of swap space to run your program.

Normally a swap partition is created during the installation of Linux.  If your swap partition is too small, you may temporarily add a swap file to enlarge it.  It's not difficult.  Let me know if you need help with doing that.  Perhaps your Mac machine has a large swap partition and that's why things work smoothly there.

@Carl Love Your post came just in the right time. I was looking for a way for doing just that.  Vote up!

In an attempt to understand what you are asking, I searched for "Bart Simpson, epicycles" on Google and found a YouTube Video which I think is related to your question.  Let me summarize.

We have a line drawing in the complex plane.  We parametrize the drawing as F(t), where t goes from 0 to 2*Pi.  Then we expand F(t) into the complex Fourier series, and interpret the partial sums of the series as depicting a planet, a moon, a moon of that moon, and so on.  Thus, the Fourier series approximation of F(t) may be viewed as the orbit of the last moon in the sequence of the moons.

In the attached Maple worksheet I show how to animate the moons as they trace the curve corresponding to F(t), which, for the sake of illustration I have taken to be a square.  The animation below corresponds to picking 7 terms of the Fourier series, and consequently the target square (drawn in green) is rendered with rounded corners (rendered in red).

Download worksheet: mw.mw

You may replace the square with any other drawing provided that you can obtain a parametric representation for it.  I understand that you are looking for a method of extracting the data for a given curve from an image, and then parametrizing it.  That's a rather complex operation.  I haven't given it much thought.  The video quoted above suggests using a traveling salesman algorithm to connect points extracted from a graph.

Added Dec 5:

Here I have replaced the square with an F-shape:

@ecterrab I updated to the latest version of the Physics package and the trouble went away.  I have added a postscript to my earlier post indicating that pdsolve() was not at fault.

@ecterrab Hi Edgardo, I have narrowed down the problem to the subscripted arbitraryfunctions business.

Specifically, in my mapleinit file I have:
`pdsolve/arbitraryfunctions` := "subscripted":

If I remove that line, then 
pdsolve({pde1, pde2, ic1, ic2});
works just fine.  Otherwise I get the error that I had noted.

The attached worksheet was executed with the troublesome line commented out in mapleinit and it shows that pdsolve works in that case.  Further down in the worksheet I insert
`pdsolve/arbitraryfunctions` := "subscripted":
after which we see that pdsolve fails.

Download mw2.mw