I have two problems really, the first problem is to find out whether the non-converging trajectories result from the true dynamics of the system (I have a local existence result, but no global existence, so maybe after all the divergence is what actually happens), and --- assuming that the divergence is caused by the lack of accuracy of the simulations (most likely my inability to use dsolve, digits, and all that) --- the second problem is to get those tricky trajectories to converge.

To investigate the behaviour near the points, I simulated a bunch of trajectories (500 to be exact) in the neighborhood of the points of interest.

This is what I get near a point in the blue region, where simulations converge:

This is what I get near a point in the red region, where simulations diverge:

Out of 500 initial conditions in a neighborhood of this "stiff" point, only a small number converged (or appear to have done).  I used the default dsolve method (rkf45) and default values for relerr, abserr, etc..

Perhaps I need a scheme to identify such initial conditions and take it from there. I'm still thinking... Suggestions always welcome.

EDIT. I should add that in both cases trajectories from "below" converge well, it is trajectories from "above" that cause problems. These are the trajectories characterized by a sharp change in direction very near the stationary points. The sharp change in direction occurs closer to the stationary point, the deeper you go into the "red-colored" region, and that is precisely the region where dsolve doesn't get convergence. You can see this from the picture I posted earlier (http://www.mapleprimes.com/view.aspx?sf=124914/419404/urq50a.jpg).

EDIT 2. Further analysis of the system revealed a condition on the parameters necessary for stability. I discuss it in a comment below.

Does anyone know if Peter Stone's procedures can be used to analyze systems of 3 (or more) non-linear ODEs ?

I couldn't find an email to contact Professor Stone --- the email given in the worksheets expired a few years ago. His webpage is still very much active since he's posted pictures of very recent hikes (nice pics by the way).

http://www.peterstone.name/
http://www.peterstone.name/Maplepgs/Maple/procdrs/DEsol.zip

I downloaded

http://www.mapleprimes.com/view.aspx?sf=119397/381650/Minecraft-Exports.mw

and did as instructed: "Click on the code button below to load the export code".

but got the following message:

3d Plots -> MinecraftError, (in unknown) external lookup of gzopen: The specified procedure could not be found.

I do have Maple 15. Any suggestions? Thanks.

@pagan 

Thanks pagan for your feedback. I should have said something about the warnings.

Do you have any suggestions in this respect?

The warnings related to the "events" are exactly what I expect. I set them up to prevent Maple from integrating beyond the region of interest.

The non-converging trajectories (red and red-ish colours) have a different kind of warning.

1. If I comment out the relerr, abserr, initstep, maxstep and use rkf45 or rosenbrock, for a quick and dirty simulation, I get something like:
   
   Warning, cannot evaluate the solution further left of -14.780539, probably a singularity
   
   
2. If I set relerr, abserr, initstep, maxstep to attemp a finer simulation, with initial points closer to the stationary points, then I sometimes get a complaint about maxfun being exceeded. Then if I raise maxfun, this message is usually eliminated, but instead I get the warning about "probably a singularity"

The maxfun warning can be eliminated by setting maxfun well in excess of the default value, but this doesn't help dsolve converge, in my experience.

The animated gif of my simulations is too large to be hosted with mapleprimes, so I uploaded it elsewhere and placed a link to it above. The preview is pretty bad, I don't know if this is caused by my browser or what. What you can see is a sharp change in the slope of q (relative to u as well as relative to time t). The red line and some of the lines right next to it are simulations that didn't converge. Any suggestions welcome! Thanks.

This is great help acer, I'd been wanting to do this for quite some time. Wouldn't life be easier if one could simply append these "play=true" options to an animation? It doesn't seem to be possible to do e.g.

plots:-display(mypreviouslycreatedanimatedplot, 'play'=true, 'continuous'=true);

I hope it will be possible to do just that in future versions of Maple, it's so intuitive.

Since I came to this post not knowing a thing about "embedded components", let me sum up here, for easy reference, what I have learned by following up on acer's link.

This is how to create the plot component. I had never done that before, indeed I had never used the palettes before.

1. Toggle the palette on the left-hand side of the worksheet.
   
2. Under the "components" tab, select the picture of a plot. This inserts a plot into the worksheet. The plot is blank and named "Embedded Plot Window."
   
3. Right-click on the plot and select its "component properties" tab. The "component properties" opens up a window named "Plot Properties" (weird).
   
4. In "Plot Expression" enter the name of your animated plot, which in acer's example is simply anim. The blank box will soon fill with the animation.

I personally detest having to use the mouse and pull down menus and stuff. The way I would set it up then is to add

   SetProperty(plotcomp,'continuous',"true");

to acer's RunMe procedure. This way it will play for as long as you care to watch. You won't need to keep pressing RunMe. And then to define the following simple StopMe procedure.

StopMe:=proc(plotcomp, checkcomp)
uses DocumentTools;
      SetProperty(plotcomp,'play',"false",'refresh'=true);
end proc:

This way, you can have:

RunMe('Plot0','CheckBox0');  # runs it programmatically

StopMe('Plot0','CheckBox0'); # stops it programmatically

Upon the first press of "enter" on the keyboard, the animation will run, upon the second press of "enter" it will stop.

Thanks again acer, great tips.

This is great help acer, I'd been wanting to do this for quite some time. Wouldn't life be easier if one could simply append these "play=true" options to an animation? It doesn't seem to be possible to do e.g.

plots:-display(mypreviouslycreatedanimatedplot, 'play'=true, 'continuous'=true);

I hope it will be possible to do just that in future versions of Maple, it's so intuitive.

Since I came to this post not knowing a thing about "embedded components", let me sum up here, for easy reference, what I have learned by following up on acer's link.

This is how to create the plot component. I had never done that before, indeed I had never used the palettes before.

1. Toggle the palette on the left-hand side of the worksheet.
   
2. Under the "components" tab, select the picture of a plot. This inserts a plot into the worksheet. The plot is blank and named "Embedded Plot Window."
   
3. Right-click on the plot and select its "component properties" tab. The "component properties" opens up a window named "Plot Properties" (weird).
   
4. In "Plot Expression" enter the name of your animated plot, which in acer's example is simply anim. The blank box will soon fill with the animation.

I personally detest having to use the mouse and pull down menus and stuff. The way I would set it up then is to add

   SetProperty(plotcomp,'continuous',"true");

to acer's RunMe procedure. This way it will play for as long as you care to watch. You won't need to keep pressing RunMe. And then to define the following simple StopMe procedure.

StopMe:=proc(plotcomp, checkcomp)
uses DocumentTools;
      SetProperty(plotcomp,'play',"false",'refresh'=true);
end proc:

This way, you can have:

RunMe('Plot0','CheckBox0');  # runs it programmatically

StopMe('Plot0','CheckBox0'); # stops it programmatically

Upon the first press of "enter" on the keyboard, the animation will run, upon the second press of "enter" it will stop.

Thanks again acer, great tips.

@Alejandro Jakubi 

and hopefully with Maple 16!!

@Alejandro Jakubi 

interesting, the maplet driver indeed does export all 3 labels. It also produces a prettier plot. Comparing like with like I exported the same plot as gif (the maplet has only a limited number of export options, with gif the prettiest, I think, for 3d plots):

with the maplet (navigate the menu):

with inline plot (right-clicking)

The good news:  the maplet driver (under standard gui) gives me better results than simply right-clicking on the plot to export. I will turn to this when I need the better-looking plots.

The bad news: involves the mouse (or a clever keyboard move) and therefore is a non-automated labor-intensive activity.

To sum up, both mouse-based methods (maplet and inline) properly display all 3 labels (START EDIT: err, no, in this case the inline right-clicking method ate part of one of the labels, this being caused apparently by over-zealous cropping. EDIT END); the command-line plotsetup(gif) "eats" the third label.

@Alejandro Jakubi 

no I haven't tried the classic GUI, which I was using until recently, because parts of my code no longer works with it, so I suspect the same problem will exist with Maple 12, though I haven't tried because I do not have access to older versions here,

thanks Alejandro,

Patrick.

The mapleprimes file uploader cannot handle long file names, if you try to upload different files that have the same first string of letters like myveryveryverylongfilename_1 and myveryveryverylongfilename_2, it will not upload the second one. It appears to have a limit of 20 characters it can handle in the name.

And even stranger is that it has posted the same image above despite getting it right during the preview!

So let me do it again. Export by mouseclicking adn you get:

@Markiyan Hirnyk 

tested epostma's code, it worked here (Windows Vista, 32bits, Maple 15, standard)

@Markiyan Hirnyk 

tested epostma's code, it worked here (Windows Vista, 32bits, Maple 15, standard)

thanks Robert for these explanations,

Christopher, I do not have the equations at hand (I'm currently traveling), but a similar phenomenon can be seen with:

plot(x^(1/2)
, x = 0 .. 1
, 'color' = blue
, 'style' = point
, 'numpoints' = 50
, 'symbol' = solidcircle
, 'symbolsize' = 8
, 'adaptive' = false
):
plots:-display(%);

plot(x^(1/2)
, x = 0 .. 1
, 'color' = blue
, 'linestyle' = dot
, 'thickness' = 3
):
plots:-display(%);

The first plot looks better with 'adaptive' set to false, but it does not "replicate" the second plot. The 'adaptive' option has no effect on the second plot.

@Robert Israel 

Great Robert, thanks a lot (I thought there was something a little odd for a small number of points),

works great,

Patrick.