Higher order numerical PDE schemes near boundaries, implementation in MATLAB Followup to my previous question. The first order scheme proved unstable for my pde: 
$$f_t + A y  f_x - B x  f_y =0$$
So I'm looking to implement a higher order scheme (using these tables). I was thinking 6th order in time and 4th in space (I'm guessing what the table gives as "Accuracy" is order?). But I don't know how to deal with the boundary conditions? I'm solving the pde on a grid with zero on the boundaries (to model the distribution vanishing in infinity), so what do I do with terms like $f(x-3h,y,t)$ at points where $x=h,2h$? Some sort of extrapolation, or another method altogether, as suggested by David in my previous question?
Here is the current MATLAB code (big thanks to David once again!)
%// Liouville equation
clear;

%// Equation Parameters:
Xmin = -10.0; %// Minimum X
Ymin = -10.0; %// Minimum Y
Xmax = 10.0;  %// Maximum X
Ymax = 10.0;  %// Maximum Y
Tmin=0;       %// Minimum time
Tmax = 2.5;   %// Maximum time
A = 1.0;      %// A parameter
B = 2.0;      %// B parameter

%// Simulation parameters:
Nt = 5000; %// Number of time steps
Nx = 150;  %// Number of X space steps
Ny = 150;  %// Number of X space steps

%// Create time and space vectors with exactly Nt, Nx, Ny elements
t=linspace(Tmin,Tmax,Nt);
x=linspace(Xmin,Xmax,Nx);
y=linspace(Ymin,Ymax,Ny);

%// Make a grid of x and y coordinates    
[X,Y] = ndgrid(x,y);

%// Calculate dt, dx and dy
dt=t(2)-t(1);
dx=x(2)-x(1);
dy=y(2)-y(1);

dtdx = dt/dx;   %// For simplicity
dtdy = dt/dy;

u=zeros(Nx,Ny,Nt); %// Initialise u matrix

%// Initial condition
u(:,:,1)=reshape(mvnpdf([X(:),Y(:)],[1.0 ,0.0],[0.25,0.25]),Nx,Ny);

%// Time stepping algorithm
i=2:Nx-1; %// Interior spatial points
j=2:Ny-1;
for k=2:Nt
    u(i,j,k) = u(i,j,k-1) - 0.5*A*dtdx*Y(i,j).*(u(i+1,j,k-1)-u(i-1,j,k-1)) ...
        + 0.5*B*dtdy*X(i,j).*(u(i,j+1,k-1)-u(i,j-1,k-1));
end

%// Animate the solution
filename = 'testnew51.gif';
for m=1:25:Nt
   contourf(x,y,u(:,:,m))
   title(num2str(t(m)))
   drawnow()
   frame = getframe(1)
   im = frame2im(frame);
   [imind,cm] = rgb2ind(im,256);
         if m == 1;

          imwrite(imind,cm,filename,'gif', 'Loopcount',inf);

      else

          imwrite(imind,cm,filename,'gif','WriteMode','append');

      end

   pause(0.05)

end

 A: The complication here, as you've pointed out, is that for a higher-order method you include more points in the derivative calculation, but near the boundary you only have one or two points.
The space and time stepping is $O(h)$ in your current method, which means that decreasing the step size by half gives twice as much accuracy, broadly speaking.
For time-stepping, you need to use an explicit method, but you can initialise it with a lower order method. For example, to approximate $f_t$, you can use the second order scheme
$$f_t\Big|^k_{ij}=\frac{f^{k-2}_{ij}-4f^{k-1}_{ij}+3f^k_{ij}}{\Delta t} $$
which is fine for $k\geq2$, but for $k=1$, we'd use
$$f_t\Big|^1_{ij}=\frac{f^{1}_{ij}-f^{0}_{ij}}{\Delta t}. $$
So using the first order order method we could calculate $f^1_{ij}$ for $f^0_{ij}$ using the first order method, and then we can start using the second order method.
For the spatial stepping it's possible to do the same thing, but because you have boundary conditions on each side of the domain, use a second order central difference method, such as (off the top of my head---check this)
$$f_x\Big|^k_{ij}=\frac{-f^k_{i+2,j}+8f^k_{i+1,j}-8f^k_{i-1,j}+f^k_{i-2,j}}{\Delta t}$$
but of course you would have to use a first order method to calculate the $f$ values adjacent to the boundary.
Ultimately this may not fix the stability problems, but hopefully it will give you more accurate solutions and will be more stable. Using a first-order method for the points adjacent to the boundary does not change the global error bounds, which is a very handy thing.
An implicit method (Crank-Nicolson, for example) will probably be more stable, but is more difficult to implement.
