@cskoog 

Chris, Thanks. It works. I have figured out the issue. Maple 6 used to have method =mgear in dsolve/numeric. This was deprecated. It is sad. 

method=mgear provides for numerical solution based on Gear's method (BDF) with numerical Jacobian. Lsode cannot be made to run in this mode (hope I am not wrong). It will be great if this can be brought back in future versions.

Just changing mgear to other stiff methods was crashing the codes because of RAM issues.

showstat shows that the procedure exists in Maple 7 (not Maple 8), but it is disabled when calling from dsolve/numeric. Perhaps someone more knowledgable (Allan Wittkopf) can comment on how to call mgear method in Maple 7 or by storing the procedure as text file and calling it in later versions.

My tests suggest that mgear was one of the best algorithms, (better than Rosenbrock type algorithms) for largescale problems (RAM issues).

@Christopher2222 

I asked them only on Friday, so perhaps they will get back to me next week. What is clear is that Maplesoft doesn't have downloadable links anymore for Maple 6.

@taro and Cskoog,

I am able to open the files, but not run the same. It may be possible that some of the commands and procedures might have been deprecated from Maple 7 onwards, not sure.

It is a big set of procedures with intentionally ignoring the warnings in multiple places, so hard to debug.

@John Fredsted 

Chinthan is in my group. I think seq command is missing

try something like 

restart;with(CodeGeneration):

> f:=[x[1],x[2],x[3]-x[2]^2];
>
> C([seq(g[i] = f[i],i=1..3)]);

f := [x[1], x[2], x[3] - x[2] ^2]

g[0] = x[0];
g[1] = x[1];
g[2] = x[2] - pow(x[1], 0.2e1);

optimization seems to find it. I think trig form may be easier to prove, I will keep trying. The equality constraint is the equation of a circle with radiaus sqrt(2). Right now, the plots make sense.

Download Inequalityproof.mws

It may be possible to combine eq2 to a form 2*sin(a*theta+b+..)^c or something like that.

I gave up on copying and pasting large amount of codes/procedures in mw format. Try mws format, it should be 100x faster for copy/paste/edits.

@ Bendesarts

Note that long time solution of oscillatory equations is hard, in particular if it has both stiff and oscillatory parts.

Also, I am a little surprised that MATLAB works for this. Can you indicate which solver you use in MATLAB? As far as I know, MATLAB does not have an inbuilt symplectic integrator, I may be wrong.

Also, thanks for this example and post. One of the things I try to teach students is to come up with simple schemes when standard blackbox solvers fail. This is a good example.

@tomleslie 

Tomlesllie, Maple does not have a direct DAE solver, but Maplesim has. Also, Maplesim can solve DAEs without exact IC for algebraic variables, Maple cannot. Also, according to my testing, Maplesim uses different approach for order reduction and symbolic preprocessing (in particular for DAEs) and I have seen it to be more robust (but more painful to setup compared to classic worksheet Maple).

@Carl Love 

Thanks for your solution. I find fsolve unreliable and a pain to work with. Fortunately Roots works.

Disclaimer: My main interest is in the roots from 0 to 1. Gauss, Lobatto and Radau roots have direct links with collocation, implicit runge kutta methods for solving stiff ODEs, DAEs and for optimal control.  Specific formula for weights can be obtained, but are not included here. My main interest is in the roots (in ascending order). Later on, I will write a detailed post on the use of these roots when time permits.

Download LobattoandRadauQuadrature.mws

@acer 

Is there a way to find Radau and Lobatto nodes and weights similarly?

Thanks

@Carl Love 

Dr. Doug Meade at the University of South Carolina developed a package in Maple (release 4 or 5) to solve BVPs using shooting method (for regular BVPs).

For linear problems (as in this case), this works good. For nonlinear or difficult problems, this may not work. Of course using compile=true might help sovler your problem faster. In addition, problems requiring solution for different parameters will benefit from this approach as well if the singularity is removed with possible transformations or taking x= 1e-6 instead of 0.

To be accurate, it is wrong to claim rkf45 to be more robust compared to BVP solvers in Maple. I can post many easy examples. Take a 2D laplace equation =0. Apply finite difference in x (keep ode form in y). Give some boundary conditions (try to guess derivatives). IVP method will fail very fast (as N increases). However, I agree that for a second order BVP with two known conditions at one boundary point, IVP methods are better most of the times.

@Preben Alsholm 

Preben I am assuming that you are doing yours in evalhf mode and are comparing it with Maple's dsolve numeric (without compile=true which runs in hardware floats). If you are using a procedural form to define delay variables, then efficiency will depend a lot on the storage mechanism used for delay variables/history.  If evalhf mode is not used in your code, then setting Digits:=16 in both the codes will enable a fair comparison. Maple may be taking 30,000 steps, not sure.

@Markiyan Hirnyk 

Mathematica cannot solve example 2. Please prove me wrong if possible. Mathematica cannot handle state dependent delays. 

From Mathematica's website http://reference.wolfram.com/language/tutorial/NDSolveDelayDifferentialEquations.html

Currently, the implementation for DDEs in NDSolve only supports constant delays.

I am trying to post the supporting text file, but don't know how. the upload button won't take it.

@acer 

Acer this is great (defining V locally for evalhf mode). Is there a similar escape route for Compiler:-Compile? I ask this because I wrote some codes in which I observe gain in Compiler:-Compile (not inlined, just multiple compiled procedures), compared to evalhf at lower dimensions. When the problem size increases, for some reason, the linearsolver (compiled mode)/one of the compiled procedures bombs out? I am worried that perhaps I am using too many arrays (Externally) in compiled code.