I agree with Wizard's comment that Mathematica isn't the best tool when computation speed and memory efficiency are very important and that you will very likely end up to use specialized code for such a large system.
On the other hand I think NDSolve
used as a pure solver is quite impressive and might potentially even be able to solve such a large system on a large enough machine (=having enough memory). NDSolve
certainly might be a good choice to learn which methods and features are important for your class of problems if you manage to produce scaled versions of your problem which show the expected behavior of the full system. Such experiments could be a good guide when looking for the right specialized code for the "real" problem if you learn that NDSolve
really can't handle it.
The trouble with large systems of ODEs in Mathematica is usually not really the solver (NDSolve
) itself but building the system of equations efficiently and handling the large output. Both can be overcome: in your case the generation of the system of ODEs in matrix form seems relatively straightforward, most probably you'll need to look into packed arrays and/or sparse arrays to create an efficient version of H
. To avoid NDSolve
to generate an interpolating function for each entry of the result vector there are several possibilities and what to do depends a lot on what you really need. You can e.g. only return the final values of these functions or use NDSolve`ProcessEquations
and friends to run the solver without creating the full output and only extract what you need. You can find examples for all these techniques in the exhaustive documentation for NDSolve
and also on this site.
As for parallel execution: many methods of NDSolve
internally incorporate the solving of linear systems of equations (and potentially other similar algorithms) which is to some extent automatically parallelized, other than that there are AFAIK no possibilities of making use of your 6 processors with NDSolve
. The best you probably could do is to run several NDSolves
in parallel in case you will need solutions for varying H
or initial/boundary conditions. If you look for another tool you will find that solvers with parallelization at various levels exist, but parallelism for differential equation solvers isn't a feature that can be taken for granted. Whatever you'll find, you should also note that a speedup of 6 with 6 processors is only the theoretical upper limit, it is usually very hard to even come close to that for real applications...
NDSolve
for large systems of ODEs. There are a lot of tricks you can use. $\endgroup$NDSolve
systems about half that large in under a minute. $\endgroup$