Solving a complex second order PDEs coupling with ODEs

Here I just discretize the dE/dz item by du[i, 9] = (u[i, 7] - u[i-1, 7]) / dx and du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx. There are no spatial terms elsewhere in the equations (u[1] to u[6]).

Are you sure you got the winding direction right? Upwind scheme - Wikipedia

I check the wiki, It looks like right :sweat_smile:

du[i, 9] = (u[i, 7] -  u[i-1, 7]) / dx # dE+/dz
du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx #-dE-/dz

Updating!
I have found a method to convert the dispersion item to the second order spatial difference item. By using central difference, the problem can be solved now(Tsit5 is faster than the ROCK4 by ~20% in my problem).

Thx for your help :grinning_face_with_smiling_eyes: