Suppose I want to numerically solve for $f(x,k)$ the one-dimensional Boltzmann equation for electrons in steady-state condition, given as: \begin{equation} \left( \dfrac{\hbar k}{m} \right) \dfrac{\partial f}{\partial x} + E \dfrac{\partial f}{\partial k} = \dfrac{f-f_0}{\tau} \end{equation} where all other symbols represent constant values.
In order to discretize the equation, I have defined the $(x,k)$ grid and I identify the coordinates of real space and momentum space with $x[i]$ and $k[j]$, respectively.
\begin{equation} \left( \dfrac{\hbar k_j}{2m} \right) \Big[f(i+1,j) - f(i-1,j)\Big] \Delta k + \dfrac{E}{2} \Big[ f(i,j+1) - f(i,j-1) \Big]\Delta x = \dfrac{f(i,j) - f_0(i,j)}{\tau} \Delta x \Delta k \end{equation} Now, the above matrix equation gives a left-hand side matrix which is not stable for arbitrary values of constants and results in very ugly oscillations.
I have seen that many people use a so-called "staggered grid" (consisting of direct/adjoint grid nodes), or use various splitting methods (even/odd or positive/negative splitting) in order to solve this problem. The more I read on these stabilization schemes, the more confused I get about the main idea. Can somebody help me with an easy-to-understand explanation or reference on the basic idea and the details of its implementation in such cases?