@Carl Love 

My computers have 64GB RAM, but that doesn't help. The issue is if you have 200 lines of codes and if you want to copy paste from one code to another or copy multiple lines of codes or edit, my mouse won't even work for 5 minutes in mw version. In mws versions, it is a breeze. I can copy 100s of lines with no trouble. Initialy my students try out the mw format, but eventually move to mws as they can't do efficient code writing/editing.

The only advantage I have found using mw is this. When you run executable codes from *.mws, the codes run and output is stored in text files (lets say). In mw format when you run this, it will open a nice window (text file) and you can see the results. It makes it easy to debug linking/compiling or other commands. I basically use mw format to run exes.

Regarding plots, that is more on the asthetics, I will pass on that. For any publications or presentations, I export data to excel and plot there.

When Maple stops releasing *.mws, I will stop paying for newer versions.

@yellowcanary 

Stick to 1D input and preferably *.mws format. You will save a lot of time and RAM when you edit and run large codes. I can never work with mw format which is slow (Java based). Maplesoft unfortunatley moved to *.mw for a dumbed down version and for doing secondary school and high school math. (Also web friendly output, platform independence, etc).

If I can be proved wrong by someone who can post the same code in 2 different formats (*.mws and *.mw) and convince me that mw format is better, I am open to change my mind.

@ I think your model is linear, so you may be able to solve a larger system of equations.

As mentioned before, lsode, gear, etc comes with limited functionality. range, parametric, compile form or events are not possible with lsode.

Did you read all the help files in dsolve and all the posts about the same?. Recently Carl gave examples for parametric output.

If you can give an example of 4 ODEs in which you use only 2 and explain what you want then you might get better answer. No one has the time to re do codes. Post the code and then ask for help. Don't post 100000 ODEs, just post 4 ODE case.

Also, output = array (datapoints) uses less RAM, so you can get output at certain points in time, interpolant = false can save memory too.

I don't want to disappoint you, but Maple's dsolve is not capable of solving large systems of ODEs as it does not use sparse methods to store or solve.  This is relevant only if use implicit=true, otherwise, it is just marching in time, so there is no linearsolve other than just marching in time as done for most solvers in Maple's dsolve.  Also stiff or any methods using Jacobian will be slow in Maple as it uses and stores a full matrix for the Jacobian instead of sparse methods.

You may be better off going to FORTRAN or C for that.

The highest I have solved are probably in the range of 500 to 1000 ODEs.

Of course if you write your own algorithm you can solve 10000 or more ODEs. 

I think you like using the GAMMA function, I understand that. If you can avoid that type of special function, C or Fortran codes will help or you have to define your own subroutines to find that.

and stiff = true? that works for this prolbem and will give you a procedure output with options to use parameters. You can use initial x value as an unknown parameter

@ The 64bit Maple2015 comes with its own compiler. Unfortunately we can use only Digits =15 when you call dsolve with compile=true. I think this needs attention from Maplesoft folks. As of now, no gain in 64bit compared to 32 bit compiler. I was hoping that 64 bit with compile will allow for higher number of digits to be used for dsolve with compile=true, may be I am missing some information/theory.

For lsode, even lsode[backfull] works for this problem. I would suggest that you use this if this problem is stiff. Backward solvers are better for stiff problems.

Also, method = mebdfi works. But won't give result at x = 0.1, but works from 0.1 to 1. 

Both of these methods are slow and won't take compile=true.

Stiff =true works if use Digits:=20. This means that you can't use compile=true option for the default solvers as well. This might have changed in recent 64 bit versions, not sure.

Your problem seems to very badly scaled. The Gamma term seems to be nearly zero till x = 0.8 and then increases. Rescaling to a different x variable might help.

@cr-ms 

Try to change the DAE solver used. Some solvers might consume too much RAM. I am not sure 64GB will help for any general case. If you keep increasing the size, you might face the same problem.

To add to this, Maple has provided lsode, gear and other methods just for educational purposes. They are not implemented in an efficient manner. (They can only solve ODEs while the original Gear or LSODE methods can do DAEs).

For stiff problems, using stiff=true or rosenbrock methods provide the same answers. For some problems, using stiff=true, implicit =true provides much more robust solver (the most robust solver in Maple based on my experience).

I would strongly suggest scaling the variables and using stiff solvers (use stiff=true). Also, try different tolerances (1e-20 or so if you are interested in very small numbers).

It is not a good idea to integrate past or near singualrity. Perhaps DAE formulation combined with events will help. If you post the mws fomrat code, I will try.

@ 

Carl Love posted similar approach. Sorry for the duplication.

Unfortunately help files are not good for this. The examples given are very easy and simple. Dealing with matrices, vectors or arrays and converting external C file to Maple is not well explained.

While Compiler Compile might work, this is not useful for future sessions and it is not trivial to use it for future purposes. 

Providing an example for converting external Fortran or C files for Newton Raphson method (not single variable, but multiple variables), Runge Kutta methods, etc will help a lot. My strong belief is that instead of adding more or more routines in Maple, Maplesoft should rather make it easy to use existing C or Fortran solvers or routines by providng more examples on wrappers, linking multiple DLL files or C files, etc. 

Sorry for the rant, of course if you were to post your Maple or C code, I can find out if I can do it. So far, I have been successful in accessing very few C files in dll form (exe format is much easier to do).

@Allan Wittkopf 

Thanks for your attention and help. Your solution though unelegant, is great as it combines fsolve in a procedural form. (In fact I was surprised that even compile = true works).  I think this is the approach you will need for DAEs as well where in you can find initial conditions for algebraic variables.

Though fsolve works for this case, it won't work for largescale problems as the ODE solution might need only 1e-6 or crude tolerances but fsolve needs more.

If I might I will add that the solvers convert the system to dy/dt = f form, not just linear (say Mdy/dt = f). If M is allowed to a function of y or at least constant, we can solve more problems.

@joha 

I would appreciate if you can provide meaningful ranges for the parameters of h. There is way to solve this without solving for dh/dt at t = 0. I would like to check my codes. If you have a reference, please mention the same.

Thanks