function[Qx_b,Qy_b]=so_diff_backward(Q,J,K,jmax,kmax,dx,dy)

%Indexing
J1=[2:jmax-1];         K1=[2:kmax-1];
J2=[3:jmax];           K2=[3:kmax];

Qx_b=zeros(size(Q));
Qy_b=zeros(size(Q));

%Difference Interior Points:                2nd order scheme
Qx_b(J2,K1,:)=.5*(Q(J2-2,K1,:)-4*Q(J2-1,K1,:)+3*Q(J2,K1,:))/dx;
Qy_b(J1,K2,:)=.5*(Q(J1,K2-2,:)-4*Q(J1,K2-1,:)+3*Q(J1,K2,:))/dy;

%Difference left and bottom boundaries:     1st order scheme
Qx_b(2,K1,:)=(Q(2,K1,:)-Q(1,K1,:))/dx;
Qy_b(J1,2,:)=(Q(J1,2,:)-Q(J1,2,:))/dy;