Numerical integration of a 2D advection equation with spacially varying coefficients. I am interested in a simple algorithm for solving this PDE for $f(x,y,t)$:
$$f_t + A y  f_x - B x  f_y =0$$
With $A,B>0$. The initial condition is some arbitrary smooth function (probably gonna use a gaussian) $$f(x,y,0)=g(x,y)$$
And $f$ vanishes at infinity.
I have very little experience in numerical methods for PDEs. I am interested in learning a simple (as simple as possible without being hopelessly inaccurate!) algorithm on how to deal with such a problem. I imagine I will be discretizing both time and the spacial variables and then replacing the derivatives with their approximations. 
The Lax-Wendorff method (http://en.wikipedia.org/wiki/Lax%E2%80%93Wendroff_method) has been suggested to me, but I can't find any entry-level materials on this, and don't even see how to apply it here (i.e. how to transform what I have into the general form). 
From what I've managed to find out, I will need to form some sort of grid for time and position. But first I would probably want to introduce dimensionless variables to replace $x,y$. Then, I will need to choose an appropriate step in the space variables and time (how do I know what to pick, so that my solution is stable?). And finally, I will be needing some sort of time step scheme, in the form
$$f^{t+1}_{i,j} = F(f^t_{i,j} \ ; \ f^t_{i+1,j} \ ; \ ... \ ; \ f^t_{i-1,j-1})$$
So that the value of the function at some node of the grid depends on its (and all 8 adjacent nodes') values in the previous time step.
In principle, this seems to be to be it, now I just need to find the appropriate discretization and scheme?
 A: The simplest time-stepping method is a FTCS scheme (forward time, centred space).
I will denote $f(x_i,y_j,t_k)$ by $f_{i,j}^k$.
This method is explicit in time, and approximates the time derivative by
$$\left.\frac{\partial f}{\partial t}\right|_{x_i,y_j}^{t_k}=\frac{f_{i,j}^{k+1}-f_{i,j}^{k}}{\Delta t}.$$
The spatial derivatives are approximated by
$$\left.\frac{\partial f}{\partial x}\right|_{x_i,y_j}^{t_k}=\frac{f_{i+1,j}^{k}-f_{i-1,j}^{k}}{2\Delta x}$$
and
$$\left.\frac{\partial f}{\partial y}\right|_{x_i,y_j}^{t_k}=\frac{f_{i,j+1}^{k}-f_{i,j-1}^{k}}{2\Delta y}.$$
(Setting $\Delta x=\Delta y$ is probably going to be helpful)
And putting it together, you get
$$\frac{f_{i,j}^{k+1}-f_{i,j}^{k}}{\Delta t}=-Ay_j\frac{f_{i+1,j}^{k}-f_{i-1,j}^{k}}{2\Delta x}+Bx_i\frac{f_{i,j+1}^{k}-f_{i,j-1}^{k}}{2\Delta y},$$
which gives this time-stepping scheme
$$f_{i,j}^{k+1}=f_{i,j}^{k}-\frac{Ay_j\Delta t}{2\Delta x}\left(f_{i+1,j}^{k}-f_{i-1,j}^{k}\right)+\frac{Bx_i\Delta t}{2\Delta y}\left(f_{i,j+1}^{k}-f_{i,j-1}^{k}\right).$$
So you can calculate $f^{k+1}_{i,j}$ using 5 grid points from the $k$-th time step. You need boundary conditions to work out the boundary values of $f$, and initial conditions for $f^0_{i,j}$, and then at each time step just apply the scheme to each interior point.
Stability is tricky in 2D. Basically, keep your spatial steps small, and your time step s smaller, and you should be alright, but because you have $yf_x$ and $xf_y$ terms you will have more problems with stability when $x$ and $y$ are large, so bear that in mind.
I hope this helps, and it should help you implement other schemes of you need higher order convergence.
