@ 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.

@JonMcLoone 

As mentioned before by me Maple does not have a direct DAE solver (it converts DAEs to ODEs and integrates as of Maple 2015). Mathematica has IDA (Backward difference based solver) for index 1 nonlinear DAEs and can solve DAEs Maple cannot solve. 

That was the only reason I tried to use Mathematica as I had to solve DAEs.

That said, my attempts to use Mathematica kept pushing me back to Maple (classic worksheet) as I couldn't handle 100s of lines (copy paste multiple procedures) in Mathematica (too slow). I would rank the regular window (mw) comparable to Mathematica in terms of speed and ease of use (copy/paste).

In addition, setting higher infolevel helps me see the algorithm behind a particular method or command in Maple (at least for most of it, not true for financial portfolio/newer packages). I don't think one can see the algorithm behind Mathematica. Please correct me if I am wrong. (For me this is very important. In fact I have learned more mathematial methods/algorithms from just browsing through Maple than reading text/papers. So, I have lost interest in Mathematica).

Mathematica used to provide access to KNITRO (not sure if it is true anymore) for optimization which got my attention for a while (KNITRO, SNOPT and IPOPT are some of the best approaches for solving large scale optimization problems. TOMLAB is a good start for those who want to see what is out there). The internal and external optimization package called by Maple is not the best and you will hit RAM limitations very fast (my guess is that zeros in Jacobian are stored and that consumes memory, moving to sparse storage would help).

When I teach, I am being pushed to use MATLAB by many, but I insist that students learn and program in what they like and I will use Maple for the symbolic stuff (generate equations/models) and use C/FORTRAN for performance.  MATLAB for control/optimization tool box.

@maple2015 

try  using fsolve first. If not, the past link might be useful.

http://www.mapleprimes.com/questions/204206-Towards-A-Compiled-Newton-Raphson

@Preben Alsholm 

Yes, theoretically fsolve should work. But efficiency of fsolve is pretty bad compared to what I posted. You can check the cpu time. Importantly for difficult problems and large scale problems fsolve may not work.

Also, one can specify abserr=1e-15, relerr=1e-13 to get more accurate solution for the dsolve solution. 

in a PDE form with well defined bcs and ics. Population balance models can be very generic and complicated (with moving particle size as opposed to constant particle size), x,y,z and time variables. So, you need to state the equations and then you can expect a Maple solution.

Of course there are some very simple models with hyperbolic/parabolic type equations with analytical solutions.

As far as I know, no one claim that they have a code to solve any population balance model with all the complexities (this includes FLUENT/CFD commerical software).

If you want to solve this, then search and try COLDAE (FORTRAN code).

For some problems you can differentiate the algebraic equations and use the equation as bc to solve this. If you want me to try this, post mws code or mws format code (Maple format as opposed math format). Your problem looks simple, so this approach might work

@tomleslie 

I checked both maple 18 and maple 2015 (not 2015.1) in 64 bit windows 7, dsolve numeric compile =true works. The 32 bit calls the watcom, but the 64 bit calls something else, I can post a screenshot for both 32 and 64 bit versions if that helps.

Regarding your point on the same compiler version, point well taken. It is trivial to write a compiled procedure (single one, not linking multiple procedures), but sadly dsolve numeric/compile won't recognize this.