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
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