What is the best way to solve finite difference equations Ax = b?

I am trying to solve a PDE using finite differences (FD) in 3D. This is the equation I am trying to solve:

Where U_0(z) is a known analytical function, L^+ is a big number O(10^5) and nu_t = kappa*z is a linear function of z. The right hand side in the equation is dependent on w so I need to iterate this equation but that is another story, for a given w my f_1 and f_2 are known.

As you see there are plenty of higher order derivatives in the equation. I implement this by constructing FD matrices (D_x, D_y, D_xx and so on) using the Knocker product in Matlab (also have grid stretching). I use stencils up to the fourth derivative and mixed derivatives are calculated as D_xy = D_x*D_y.

Now to my question. What is the best way to solve this equation in Matlab? I started with mldivide (x = A\b). However, I found that when the size of the problem increases the speed of mldivide start to become really time consuming. The goal is to do a calculation with around 200x80x80=1280000 grid points. I had a similar problem before where I used bicg (Matlab function) instead of mldivide and could boost the performance and still have a robust solution method. Now I though it would be more or less the same but this time bicg does not perform as I expected. It take much more time but to compute than I hope for even for a small problem (around 50x30x30). I should add that the matrix is not positive definite or symmetric and that I tried GMRES without any success. So I started to think about preconditioners and tried ilu which is successful for smaller problems (around 30x20x20 up to 50x30x30) but then it also becomes really slow. And for the small problems mldivide is outperforming bicg+ilu anyway. Then I though of physical prconditioners (where you use the actual physical problem to construct a preconditioner). But I am not sure how to do this, I tried to get ride of the "coefficients" so that the equation would be:

and then use that as an preconditioner but it then bicg takes even more time and it ends with me breaking the call. Lastly, I considered the possibility of using multigrid methods but I think that will take too much time to implement if there isn't any code that is already implemented. Anyone that has any suggestions on how to solve this equation more efficient? I have added an example A matrix and b vector for a system that is 80x50x50 in this dropbox link:

https://www.dropbox.com/s/ygqqyuidxv98y4i/FDsys.mat?dl=0