I have been transitioning my MATLAB code to Julia and managed to get working a sort of analogue to my MATLAB code.

My PDE system has 2 coupled (1+1)D nonlinear PDEs:

-\mathrm{i}\left(\frac{\partial \psi}{\partial X}\right) = -\frac{\mathrm{sgn}\left(\beta_{2S} \right)}{2} \! \left(\frac{\partial^{2} \psi}{\partial \tau^{2}} \right) + {| \psi|}^{2} \psi,\\
-\mathrm{i} \left(\frac{\partial \phi}{\partial X}\right) = -\frac{s_{4}}{2} \left(\frac{\partial^{2} \phi}{\partial \tau^{2}}\right)+{i} s_{1} \! \left(\frac{\partial \phi}{\partial \tau}\right) + s_{2} \psi \phi^{*} {\mathrm{exp}}[{\mathrm{i} \left(s_{3} + q \right) \! X}] + s_{5} {| \psi|}^{2} \phi.

In MATLAB, we used a UPPE style method to transform these equations into a set of coupled ODEs by taking the Fourier transform, solving the equations, and then inverse transforming them to find the original functions \psi and \phi.

The ODEs look like this:

-\mathrm{i} \left(\frac{\mathrm{d} \varphi}{\mathrm{d} X}\right) = \mathcal{F}\left[\frac{s_{2}}{2} \phi^{2} {\mathrm{exp}}\left({-\mathrm{i} s_{3} X}\right)+\left({| \psi |}^{2}+s_{5} {| \phi |}^{2}\right) \! \psi \right] {\mathrm{exp}}\left({-\mathrm{i} \alpha X}\right), \\
-\mathrm{i} \left(\frac{\mathrm{d} \theta}{\mathrm{d} X}\right) = \mathcal{F}\left[s_{2} \psi \phi^{*} {\mathrm{exp}}\left({\mathrm{i} s_{3} X}\right) +\left(s_{5} {| \psi|}^{2}+s_{6} {| \phi|}^{2}\right) \! \phi \right] {\mathrm{exp}}\left({-\mathrm{i} \Omega X}\right),

where \mathcal{F}\left(\psi\right)=\tilde{\psi}, \mathcal{F}\left(\phi\right)=\tilde{\phi}, \tilde{\psi}=\varphi {\mathrm e}^{i \alpha X}, \tilde{\phi}=\theta {\mathrm e}^{i \Omega X} and \Omega is related to the initial conditions and parameters of the system.

I have solved the system in Julia with some simple initial boundary conditions using a variety of different solvers in DifferentialEquations.jl (I found DP8 to be the most performant in my limited testing).

After all this, my question is am I going about this the right way in Julia? I have no experience with the PDE solvers available in Julia so I am not sure if they will be more performant than the method described above. I know this might be really hard to judge from what Iâ€™ve posted and I should probably just try both and compare, but if anyone has any insights into an approach that could prove more efficient and performant that would be great.

Thanks in advance!