@ 

While not a specialist of Maple, but simply as someone who can (in general) deal with the problems he faces, I am used to using Maple for a long time.

And over all those years I have felt some "semantic shift" in the strategy of development : more and more often "derived products" (Maple TA, MapleNet, MapleSim, Embedded components, Interactive documents and so on) are promotted at the expense of Maple itself.
I had recently questions about the 'compile=true' option of dsolve and its compatibility with the 'parameters=[...]' feature of it ... and the answer I was given was that it was the kind of things MapleSim was doing well
In the same vein, some of my colleagues recently took a 3 days training course about Maple ... but not far from 50% of the content was about those new features.


This is probably the consequence of the race between Mathematica (which was rather in advance in providing  "engineering solutions") and Maple.
Nevertheless, it would be damaging that the necessary attempt to fill this gap leads Maple to push into second place its core business and undermines trust.

I hope I will have the pleasure to read you again on other topics

@ 
I am really impressed by your analysis of lsode.


In fact I rarely use Maple for solving ODEs numerically. As I do a lot of statistical studies my "reference" tool for numerical purposes is R (https://cran.r-project.org)
Before using Maple to solve the problem I presented in my primal question, I had solved it in R (of course the equations were not written in the same form, but they represented basically the same physical system) by using the function lsode from the deSolve R-package.

In the package deSolve, it is written
The R function lsodar provides an interface to the FORTRAN ODE solver of the same name, written by Alan C. Hindmarsh and Linda R. Petzold.
(lsodar allows you to invoke lsode, lsoda, lsodar, ... and, contrary to Maple, allows use of events)

The solver to wich this sentence refers is part of the Fortran library ODEPACK
https://computation.llnl.gov/casc/odepack/odepack_home.html

I naively thought that Maple had used the same strategy, in the sense that : the function lsode provides an interface to some function of the same name ... 
For instance from  ODEPACK (Fortran) or SUNDIALS (C, C++) which are highly respected ODEs solvers libraries.
This because I thought it would have been a bit stupid to redevelop something that works well somewhere else.

But, as the head of the mail the Maple support team sent me (« The problem is that lsode is a direct translation of...), and as your sound analysis seems to suggest too, Maple has probably coded its own version of lsode  ... and maybe has done this in a way that is not above reproach (?)



Tanks again for your definitely exciting analysis.




Here is the "in extenso" version of the Maple team answer 

« The problem is that lsode is a direct translation of an involved Fortran implementation that uses common blocks (aka global variables) to store critical data, and needless to say the exact nature of these is not documented.

The issue is simple to describe:

Say you obtain a solution 'sol' to an ODE problem using LSODE, and you do this:

sol(1);

sol(2);

Then do the same thing, but then do this:

sol(2);

One might expect the two values of sol(2) to agree, but they will not, and here's why:

In the first sequence, the lsode solver is called to integrate from 0 to 1, and obtain values at 1.

For efficiency, when the value at 2 is requested, we do not integrate all the way from zero again (though we probably should) but instead use the value at t=1 as a starting point.

The problem is if we directly request the value at t=2 instead, the integration is highly unlikely to choose t=1 as a natural point to take one of the integration steps, so the computation process will be different, and the values will differ.

The solvers where it was possible to see what is actually going on (how the integrator works) were updated to store the last step before a modified step (i.e. modified by the endpoint) and resume the integration from there when another point is requested.

So as an example, in integrating from 0 to 1, we might evaluate at the following time values:

t=0, 0.11, 0.207, 0.292, 0.35, 0.38, 0.41, 0.48, ..., 0.902, 0.949, 1.00

where in effect the last two evaluation points were altered from the natural step size of the solver as driven by the error control so that the endpoint would be hit. The altered solvers save their state before the alteration, effectively at 0.902, and resume from there once a request is made for a point past t=1, so for example you might see this:

t=0.902, 0.97, 1.03, ...

The main solvers are even more intelligent in how they do this, as there is no repetition of computation as they store a antural interpolant, so for a request at t=1 they would take the steps .902, 0.97, 1.03, then simply interpolate the t=1 value over the 0.97..1.03 region.

 It was not possible to do this for lsode»

@ 

I'm happy to read a new contribution after some months of calm.

Personally I had drew a line under this question for two reasons.

First one is : despite the many contributions I read I did not find "the" answer I was looking for ; a lot of replies were very interesting but "off track" in the sense that they did not answer what I had asked.

The second is less subjective : I finally received an answer from the Maple hotline that (for me) definitely ends the debate.
As I am living in France I submitted the problem to the french support team, which did the same to Maplesoft Waterloo, which in turn provided an answer to the french support team which forwarded it to me, quoted by a remark in french.

The answer from Waterloo begins by

" The problem is that lsode is a direct translation of an involved Fortran implementation that uses common blocks (aka global variables) to store critical data, and needless to say the exact nature of these is not documented ..."

(the highlight is mine, see below)

Then there is a lengthy explanation of how lsode performs evaluations of the procedure solution for given values of the variable (needless to say that I can provide you the total answer if you are interested)

The end of the mail is in french, so I translate it for you
"Due to this problem, the lsode command will be probably deprecated in future Maple versions"


Let me know if you want to read the complete response from the support team.
My sincerest thanks for your involvement.



About the highlighted text : one can ask whether Maple acted judiciously by making avaliable a wrapper of a Fortran function, that it itself recognizes as not documented ?
This personal remark finds an echo in your last sentence :"If this is true, then it is unfortunate that Maple claims to implement standard methods, but implements the same in an incomplete manner."

@Carl Love and the owner of the line and of the Maple license. So I usually find easier to use the same email adress when I use his personnal Maple license or to access Maple Primes.
I do hope there is nothing inappropriate, don't I ?

@Carl Love Absolutly right : I do not know what I was thinking

I was so focused on the error message that I did not realize I had permuted the indices

Sorry again : if it exist some Darwin Price here I deserve it !!!

@Carl Love Don't worry about the difference between my name and "sand15" : it depens on whether I access Maple Primes from my personnal computer or from my professional network.

@Carl Love  I am not sure you have correctly read my question :
 I have written  myVal := userData[thisNode+1, 1]  
 and not           myVal := userData[thisNode+1] 

@Carl Love 

I'm interesting in the limitations of cmaple (I run Maple on a Windows PC) : what are the features that are correctly supported and what are those which are not.

For example, plot commands, whatever they are, are not "correctly" taken into account by cmaple (I run Maple on a Windows PC) in the sense that a simple command like 

plotsetup(pdf, plotoutput=`MyPlot.pdf` ...)   #the type of the output file doesn't matter
plot(x,x=0..1, color=red, [+ some extra informations for a pretty graphic] )

even if it doesn't generate any warning or error, produces a poor graphic : no options but a few of them are handled by cmaple.

Another example which I encountered recently is the following 
Here are different syntaxes for using fsolve which seems to return the desired answer (even if only the three last of them seem to match what is written in ?fsolve[details])

fsolve(sin, -1..1);
fsolve(sin(x), x=-1..1);

fsolve({sin}, {x=-1..1});  # plus (op(%)) if needed
fsolve(sin(x), x, -1..1);
fsolve([sin], [-1..1]);        # plus (op(%)) if needed

In my maple worksheet I use  fsolve([sin], [-1..1]), which returns [0].
But the corresponding mw file, when run from cmaple returns an error I was able to fix by writting 
fsolve(sin, -1..1) (which luckily works correctly in the Maple worksheet even if it is not declared as a legal syntax in the help page [if I am not mistaken])


So I suspect some features, or some "new" syntaxes, are not correctly handled by cmaple : maybe these drawbacks have been already reported and listed somewhere ?

Beyond this, I am having problems of computational order such as memory inflation or large cpu times.
A way to reduce this later could be found in using the Thread package, or even the Grid facility : but is it possible to 
use them if I run the code in a batch mode ?

These are the types of questions for which I was looking possible answers

@Doug Meade 

Thanks Doug

I thought I read somewhere in this same manual that T[1]:='T[1]' was the way to remove T[1] (which works well in fact).

But you're right I should have Read The Famous Manual more closely.

@Carl Love 
(shame on me)

thank you Carl

(once gain, it will become a habit)

Further to my previous answer : there were still other typping errors somewhere (theta is misinterpreted).

I rewrote your PDE system and try to solve it
This error message is returned

Error, (in pdsolve/numeric) initial/boundary conditions must be defined at one or two points for each independent variable

It means you have insufficient boundary conditions, in other terms your problem is found to be ill-posed.

Without more precision concerning the physical problem, and more precisely the [xi, eta] domain of definition (you write IBC = ... but do you have BC and IC or only BC ?), I am not capable to say more. 



Download PDE.mw

The problem is on your third line : xi := 0 !

@Carl Love  perfect

@Carl Love 

 
Your addition of  expand does not fix the "problem" I raised previously.
But, at least to me, it is anecdotal, I am completely satisfied with your answer

@abcd another typing error ... "too much ardour on my part"