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