I used the option output=listprocedure in the dsolve() command. This means that dsolve() will return a list of procedures for calculating all the values of interest. If you display the output of the dsolve() command you should find that it shows a list someting like

res:=[ t=proc(t) ... end proc,
          u[1](t)=proc(t) ... end proc,
          u[2](t)=proc(t) ... end proc
          u[3](t)=proc(t) ... end proc,
          u[4](t)=proc(t) ... end proc,
          v[1](t)=proc(t) ... end proc,
          v[2](t)=proc(t) ... end proc,
          v[3](t)=proc(t) ... end proc,
          v[4](t)=proc(t) ... end proc
     ]

so v[1]|(t) can be computed by accessing the rhight-hand-side of the fifith element of this list, v[2](t) by the sixth element and so on: hence v:= rhs(res[5+j]) for j=1..4 gets the 6th..9th elements of this list, corrresponding to the procedures for computing v[1](t), v[2](t), v[3](t), v[4](t) respectively.

Hence the statement

v:= rhs(res[5+j]);

just steps through the relevant entries in the list 'res' returning the procedures to compute v[1](t), v[2](t), v[3](t), v[4](t) rspectively, and assigns each of these procedures to v, just to make subsequent computations "simpler". There are many other ways I could have selected these procedures, and many other assignments I could have made  - many of whihc would be "neater" and "prettier!"- but at the time I was being I was being "quick, dirty and effective"

In order to reduce the "thickness" of the line being plotted, you have to bear in mind that my original is actually plotting "points". To reduce the apparent "size" of these points, one has to select a symbol type (asterik, box, whatever) and then st the symbolsize. If you substitue the plot stages in my original with

  tmax:=5:
  for j from 1 by 1 to 4 do
      v:= rhs(res[5+j]);
      la||j:= plot([seq(`if`(v(t)>0,[t, 5*j],NULL),t=0..tmax,0.01)],
                        style=point,
                        color=blue,
                        symbol=solidcircle,
                        symbolsize=2
                      ):
      lb||j:= plot([seq(`if`(v(t)<0,[t, 5*j],NULL),t=0..tmax,0.01)],
                        style=point,
                        color=red,
                        symbol=solidcircle,
                        symbolsize=2
                      ):
  od:
  plots[display]([p1, seq([la||j, lb||j][], j=1..4)]);

then you will get what "seem" to be *nice, thin* lines for the "gait" portion of your plots.

The only way I can think of to fake the tickmarks on the x-axis is to supply a specific list of tickmarkValues, as in changing the display(0 command in my original to

plots[display]([p1, seq([la||j, lb||j][], j=1..4)],
                  tickmarks=[default, [0=0, 5=0, 10=0, 15=0, 20=0]]
                );

It is possible to "spoof" the y-axis labels by using the 'tickmarks' option: although this is trivial to do, the aesthetics may be "unpleasing" becuase one ends up with an extenced y-axis either at the top or bottom of the resulting plot. I can think of ways to fix this generally, but it would require woring out maximum and mininum values on all the v[j](t) plots and then spacin them appropriately - and frankly this is the kind of "cosmetic" issue which really doen't interest, so I can't be bothered. I just eyeballed some appropriate values in appropriate places which make the final graph  "look good". Achieving the same effect for arbirtarty values occurring in the v[j](t) could be achieved, but would be more complicated.

The attached should fulfil all of your immediate requirements

gaitPlot2.mw

NB If the final plot does nt *seem* to have the "gaits" plotted, then stretch it in the vertical direction - if this is unacceptable then increase the 'symbolsize' value in the plot() commands

I can't help feeling that there must be a better way to do this, but if you don't mind "ugly" then the following will work

gaitPlot.mw

since I don't have access to MapleSim, I can't really test, so the following are just "observations" based on no evidence

1. MapleSim is still solving a set of differentisal equations numerically: there are only a certain number of ways in which this can be done, and I very much doubt that Maplesim contains a DE solution method which isn't available within Maple. In fact I would expect that both use exactly the same underlying code for ODE solution
2. This could still be the *visual* problem of too few points on each orbit, so I would try to determine this first - either fix the number of orbits, and increase the number of points; or plot points only to see whether your original result is an "illusion" caused by plotting the lines between the points

1. Eq2 still contains a function theta(x) - I assume this should be theta(x, xi)
2. Now that you have a system of 2 partial differential equations, you have to use pdsolve() rather than dsolve() to obtain a solution. Note that the syntax for pdsolve() is different from dsolve(): you will have to fve the equations and boundary conditions as two separate sets/lists - see the pdsolve() help page
3. You are still using a function definition (an -> operator) within your boundary conditions, but the syntax doesn't make any sense. I think when you write D(f)(xi,N)->1, what you mean is D(f)(xi,N)=1 - but of course I'm guessing, and I may be wrong. Whatever you mean, D(f)(xi,N)->1 is always syntactically wrong. Check the help page for the -> operator: is that really what you mean
4. Where you have differential boundary conditions, the form is incorrect because you have to specif, which variable you are differentiating with respect to: for example D[1](f) will differentiate with respect to the first variable in f, and D[2](f). Suggest that you read the help page for the D-operator
5. Even correcting the above, I don't think that you will have a valid set of boundary conditions because in every single one, the first argument is symbolic, and the second argument is numeric. Essentially this means that you have no boundary conditions on the first variable
6. Confusingly, the first variable in all your boundary conditions has been labelled as xi, whereas in all of your function definitions, you have the first variable as x and the second variable as xi. This isn't necessarily an error, but the consequences are confusing (at least to me), so I suggest you check these definitions very carefully

1. In the definition of your equations (eg eq1), the funtion f() appears to be a function of one variable, since it is generally written as f(x). Hwever the last term indicates that the function f(x) should be differentitated wrt x, and then wrt xi. Differntiating a function f(x) wrt some other "random" variable (ie xi) is going to produce zero. Are you sure that you want this?
2. eq2 is at least consistent in that both f() and theta() seem to be functions of one variable, however
3. in the definitions of boundary conditions, both f() and theta() seem to be functions of two variables. You really need to decide whether f(), theta() are functions of one variable or functions of two variables - until you do this, there is little I can do to help
4. In the definition of boundary conditions, you seem to *think* that you can define some kind of "limiting" behaviour - so far as I know, you can't. You need to write explicit equations (in either one or two variables)
5. I agree with Preben - change the name of your worksheet

I do not know how you have managed to make so much of your worksheet essentially unreadable - and it is too big and too slow for me to debug. However, if I painstakingly use copy/paste to get the first 10 points of your data, and look for a fit, then the problem is trivial, as in

restart:
with(Statistics):
Xdata:=[12, 15, 18, 21, 24, 27, 30, 33, 36, 39]:
Ydata:=[627.812805, 626.780884, 626.8172, 626.789246, 626.75708, 626.701111, 626.702881, 626.672974, 626.643188, 626.584778]:
LinearFit(c*t^2+b*t+a, Xdata, Ydata, t);

which returns

Now I'm pretty sure I can do the same thing with your complete data set, provided that you you make this data set available in a form which I can read. I woud suggest two simple text files, one for the xdata and one for the ydata. The individual files can be tab-separated, csv, or pretty much anything legible in a standard text editor

If I add the line

map(combine, ListTools[PartialSums](convert(U, list)))

in order to get the partial sums within your list it is obvious that the algebraic term in each such sum increases by Pi⁴/2³ or roughly 12x.

In other words if I write the n-th term as n_k*x²+cos(x), then the next term will be (roughly) 12*n_k*x²+cos(x). No way this is going to approach cos(x)

Using the DirectSearch package returned a solution for all loop index values, although

1. It complains a fair bit about complex/non-numeric values
2. I had to increase the function evaluation limit, (specifically when the loop index=8)
3. Back-substitution of the results obtained showed that the first two equations 'A' and 'B' were solved pretty well for all loop indexes, but the final two - well, not so good

Use the attached at your own risk!

fsolveProb.mw

Again I suspect that you are being confused because you are not plotting the solution points, but also the lines in-between them: see my previous post on this subject

Redo this plot with whatever option you need to replicate the 'style=point' option in Maple's plot command: do you still see a "drift"

Consider the toy example

plots[implicitplot]([x^2+y^2=1], x=-1..1, y=-1..1, numpoints=40, style=point);

Somewhere around 40-50, I am convinced that this is a circle - however if I draw the lines between the points using

plots[implicitplot]([x^2+y^2=1], x=-1..1, y=-1..1, numpoints=40);

then I no longer see a circle, I see a definite polygon. When drawing the lines between points, then my eyes need about 400 (try it) points before I see a smooth circle.

Now apply this insight to your second example above: compare the output of the commands

plots:-odeplot(res,[u[1](t),v[1](t)],0..tmax, scaling = constrained, numpoints = 10*tmax);
plots:-odeplot(res,[u[1](t),v[1](t)],0..tmax, scaling = constrained, numpoints = 10*tmax, style=point);

The points being plotted are exactly the same - but the graphs look completely different!

The toy example above illustrates that I need ~400 points on a simple closed curve in order to interpret it as a "smooth" curve.Your eyes may need more (or less). Since your example has ~25periods (or cycles) then I would expect to need ~10000 points in order to fool my eyes. This corresponds to the value of 100*tmax

From previous entries on this thread you can plot

"log ||E||_h as a function of log h"

In order to plot a "least squares approximation|" then you mean that you do not want to plot

"log ||E||_h as a function of log h"

but rather you have some arbitrary function for your dependent variable - lets us say a*log(h)+b, and you want to determine values for the parameters a and b in thi function, to obtain the best possibe fit in a least squares sense???)

Makes no sense - if you already have a set of pairs [log ||E||_h, log h] then why would you want to produce a set of pairs ["some approximation of log ||E||_h", log h].

Don't get me wrong, this can be done. However it does reuire you to determine the function "some approximation of log ||E||_h", which I used earlier - one just has to know what this function is!!!

For example if you said that you wanted to fit the function

a* sin ( log(h) ) + b*1/arctan(log(h))

then I could do it. But whihc function are you trying to fit in na least squares sense (and why????)

Like I said - you original question was ambiguous: and I gave you a solution on what you *might* mean. I don't think it is too difficult to see what was being calculated, based on your original equations.

BTW I am unconvinced by your statement that

"They are all spheres."

The last pointplot in your sequence looks nothing like a sphere: in fact it is rather like the "point" version of the first "surface" plot you produced using my interpretation???

Sometimes, in order t get "precise" answers, you really have to ask precise "questions", cos if you don't know the question, then I can guarantee that I won't know the answer ;-)

Preben is correct in that if you know that the oscillator is asymptotically stable, then it makes no sense to compute more than one period. Hold this thought through all of the following because it is important

In your original question you stated that "the plot didn't stay in the limit cycle", so I was looking for evidence of departure from the limit cycle, which I didn't find, but I didn't want to restrict myself to one orbit - just in case the sytem "popped out" of the limit cycle at some point

There was no evidence that system ever left the limit cycle, and my conclusion was that you were being confused by the fact that, as tmax increased, you were only computing three or four points on a cycle, then joining the dots, which look like a triangle (or square), and nothing like an ellipse. Hence my rather loose use of the term "subsampling", and why the apparent problem could be solved by increasing the numbers of points which are plotted. After all plotting three points on an ellipse (and joining the dots) will give you a triangle - but plot 100 points on an ellipse (and join the dots), will look fairly elliptical! So raising numpoints as tmax as increased would work: quick and dirty

However this argument does beg the question - how does one compute the "fundamental period" for such an oscillatory system, other than experimentally? If I play with tmax on your original problem, then a single period occurs somewhere between tmax=6 and tmax=7. If the oscillator is guaranteed to be asymptotically stable then it makes no sense to compute the solutions for any value of tmax greater than this: and. for a single cycle, there is no gain in raising numpoints: the default value will provide a "reasonably smooth" ellipse

What this means is that your question about how to "tune numpoints correctly" is the wrong question! The question should be, "How do I determine the fundamental period?" - and other than "experimentally", I don't know the answer: rather obviously one is looking for the case where

[x(t+tau), z(t+tau)] = [x(t), y(t)]

If one grants the status of an (arbitrarily accurate) independent variable, then this "roughtly translates" as looking for the condition

[x(t2), z(t2)+epsilon] = [x(t1), z(t1)]

for some "small" value of epsilon. In which case the fundamental period is t2-t1. However this does not really help in determining the periodic condition, since it depends on the user-selected value of epsilon.

I would actually be interested in an efficient means of determing the fundamental period of an arbitrary oscillatory system, cos right now, I can't think of a way to do it :-(

Your original request is a bit ambiguous(?). I interpreted it as

restart;
with(plots):
xv:=(p1,p2,p3)->2*p1*p3/(p1^2+p2^2+p3^2):
yv:=(p1,p2,p3)->2*p2*p3/(p1^2+p2^2+p3^2):
zv:=(p1,p2,p3)->(p3^2-(p1^2+p2^2))/(p1^2+p2^2+p3^2):
implicitplot3d( [ xv(m,n,p), yv(m,n,p), zv(m,n,p) ],
                      m=-2..2,
                      n=-2..2,
                      p=-2..2,
                      style=patchnogrid
              );

which does produce a "nice" plot

Normal procedure n this forum is that you post your code and someone here will try to fix it. Shows that you have made some kind of effort.

By the way I should have pointed in my ealier post that the the Runge-Kutta method only applies to initial value problems, and cannot (in general) be applied to boundary-value problems, which yours is.

So the results which you posted, claiming the uses of a Runge-Kutta algorithm are very suspicious!