We'll actually be releasing a set of numerical codes for the GMMNLSE for MATLAB in the next couple weeks. If you're interested, pleas send me a note or reply, and I'll send you a link when we release it. 

Pretty straightforward use of MATLAB's gpu functions leads to about 2 orders of magnitude faster simulations. We then implement a parallel algorithm that breaks up the longitudinal integration into a large number of independent computations. When these are iteratively corrected, there tends to be another OOM speed-up, especially for highly multimode sims. 

(BTW, agree with Michael's comment too, of course - especially his first suggestion).

Thank u Mr. Logan G Wright.  I am approaching to solve this with RK fourth order method but still facing to solve the problem as for high picosecond pulses programme runs very slowly. I am also facing some problem regarding coupling coefficient . I am very much excited to get your numerical code. Please please send me the link. It will be very much helpful to me. Waiting for your response .

The fourth-order RK is, in general, an overkill for this problem, which means you'd be wasting your time to operate in an UNNECESSARILY high level of accuracy, and it might cause instabilities.

The reason is that the traditional symmetric splitting is a second-order scheme. This means ANY numerical scheme that is used in a single propagation step should NOT be more than second-order --- otherwise you are wasting computer time.

You can use our code that is attached below, and also published as Wang, et al, JOSA B, 30 (11), 3064 (2013). 

On top of that, an EXTREMELY helpful technique is adaptive step sizes (see Sinkin, et al, JLT 21(1), 61 (2003)). Using adaptive step sizes --- which is built on Richardson Extrapolation --- you can not only formulate a third-order scheme, but also optimize the step sizes using the local numerical errors (which come free with Richardson Extrapolation). I myself have tried this technique, and there is simply such a huge improvement comparing with codes using constant step sizes. 

http://photonics.umbc.edu/software/HME_paper.zip

https://www.osapublishing.org/josab/abstract.cfm?uri=josab-30-11-3064

https://www.osapublishing.org/jlt/abstract.cfm?uri=jlt-21-1-61

Thank u vary much Mr. Wang. This is really nice programme with high accuracy but this is not what I exactly want. This programme is for single mode GNLSE but I want it for multimode scenario. Also this adaptive step size method can be used to make the programme for multimode. I'll try it. Thank you very much for sending this file.

https://github.com/WiseLabAEP/GMMNLSE-Solver-FINAL

There's our code. It can run on just a CPU, but you are strongly encouraged to use a GPU if you aren't ~infinitely patient.

One detail is that we have found the adaptive step size is less useful for MM propagation than for SM propagation due to the fast mode beating. You end up needing to stick to a very fine longitudinal resolution even when the characteristic nonlinear length is pretty long. I suspect there are good ways to improve on this. Since for telecoms (or lasers, if you consider many round trips), people would like to simulate many kilometers of fiber propagation, having a way to model the mode beating (and sometimes linear mode coupling) correctly while using a much coarser resolution could mean huuuuge speed ups. A new preprint along these lines was uploaded by Conforti et al. just this week https://arxiv.org/abs/1707.09223 

The algorithm we use does provide a kind of adaptive step, in terms of the number of iterations that are applied in each 'super step'. However, by and large we know there are a few additions that might provide some speed-up. Ultimately, we figured better to keep the code as simple as possible. I'd be keen to hear about any innovations, though - particularly with respect to getting something like the awesome benefits of the adaptive step size per Shaokang in MM.